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N OWERCLATURE 

DESCRIPTION 

Area 

Accumulation parameter. Eg. (34) 

Bassett unsteady memory force 

Airfoil chord length ~ 

Droplet drag coefficient 

Airfoil drag coefficient 

Airfoil lift coefficient 

Airfoil moment coefficient 

Viscous drag force 

Total.airf oil collection efficiency, 

Fq. ( 28 ) 

Froude number, Eq . (6) 

Gravitational acceleration constant 
Airfoil projected height 
Inertia parameter. Eg. (5) 

Trajectory similarity parameter, Bq . (24) 

Modified inertia parameter, Eq. (14) 

Roughness height 

Length of ice growth 

Mass of ice accretion 

Mass of water droplet 

Apparent mass force 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Pressure graclient force 

Droplet Reynolds number. Eg. <7) 

Airfoil surface radius of curvature, leading 
edoe radius 

Droplet free stream Reynolds number. 

Eg. (12) 

Effecti ve radius of curvature. Eg. (38) 
Airfoil surface arc length 
Air temperature 
T i me 

Free stream velocity 
Local flovfield vieJocity 
Cummulative volume, percent 
Hori 2 ontal and vertical coordinate 
Airfoil angle of attack 
Airfoil angle of attack vh«.n iced 
Zero lift angle of attack 
Impingement efficiency. Eg. (25) 

Approximate drag law exponent 
Droplet diameter 

Angle between surface outer normal and 
vertical 

Dimensionless droplet position 

Impingement angle 

Liquid water content of the cloud 

Ratio of droplet trajectory to Stokes law 
trajectory. Eg. (15) 


V 


Absolute air viscosity 
Air density 
Droplet density 
NondimenEional time 

A.nqle between the droplet trajectory and outer 
surface normal at impact 

Anqle between the droplet trajectory and the 
vertical -at impact 

Subscripts 

m Model 

t Total_conditions 

o Initial condition 

Superscripts 

. Derivative with respect, to nondimensional time 


y 

p 

0 

T 


Vector notation 
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I. INTRODUCTION 


The scientific study of aircraft icing began in, the 
1920 's when aircraft were first relied upon tor dependable 
transportation and national defense. Recently the 
increased utiliy of general aviation aircraft and 
helicopters has resulted in an increased potential for 
unfavorable encounters with ice. Advances in avionics has 
made instrument navigation very reliable and sufficiently 
inexpensive to enable this equipment to be within reach of 
most qeneral aviation aircraft. 

The a.e.nQcynamic penalties associated with flight into 
known icina conditions are well known; a sharp dra t.-.se 
and a reduction of maximum lift coefficient. However 
avoiding icina by remaining on the. ground when such 
conditions are predicted results in an economic penalty of 
loss of aircraft usefullnesc which is not easily accepted. 
The ph-ysical processes involved in aircraft icing, and 
therefore the solutions to the icinq problem, are very 
com plex . 

Aircraft icing occurs when an aircraft penetrates in 
flight a cloud of small super cooled water droplets. A 
portion of these droplets impinge upon the leading edges 
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of various aircraft components risultino in the qrovth or 
accretion of ice. The accretion of ice ami its effect on 
the aircraft is a very difficult problem requiring the 
expertise of many areas of science and engineering. 

However most of the pro blem falls into one of two 
categories; thermodynamics or aerodynamics. 

The thermodynamics of aircraft icing deals with the 
process by which the droplets which impinge on the surface 
change from the liquid to the solid phase. Two types of 
ioe accretions can be identified and these are depictea in 
figure 1. Rime ice forms a relatively streamline shape 
extendina into the oncoming air, while glaze ice is 
characterized by the double horned shape. Table 1 
summarizes the conditions under which each type of ice may 
be expected . 


TABLE 1 ICE F0R1ATI0N 



Rime 

Glaze 

Liquid Water Content 

Low 

High 

Air Temperature 

Low 

Near Freezing 

Plight Velocity 

Low _ 

High 

Preezino Fraction 

One 

Less Than One 

D IX) p lets Freeze 

On Impact 

Flow On Surface 

Ice Color 

White, Opaque 

Clear 

Ice Density 

< 1 qm/cc 

1 q m/cc 


Rime ice occurs at low air temperatu 
water contents (the concentration of 
free stream) and low flight speeds. 


res and at low liquid 
water droplets in the 
111 rime icing the 
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drop^Pts frppzp oTi impact and a qood approximation to this 
orosth can be made by ncjlcctin., all thermodynamic effects 
Ml . Olaxe ice occurs at temperatures slightly belo. 
freezino and ,t relatively nigh liguid .ater contents and 
hion flight velocities. An analysis of glaze ice 
accretion must include the proper thermodynamic modelling. 

a^'Todynamic wind tunnel test of a 
Simulated ice shape [2] are shovu in figure 2. Urge 
increas£^i.n__or^ ami a reduction in'naxi.um lift 
coefficient are shoxn for both types of ice. I^d 
air. oils are difficult to analyze due to the severe 
surface rouahness and lartje zones of separated flow which 
result fro. the irregular shapes of the ice accretions, 
only empirical .ethods are currently available to predict 
this dearaJatio)! in jerformanec. 

T»o npnronchcr to the aircraft icing problem are 
available. The fir.at- method is to prevent the ice from 
forming (ant-icing) or to renoye it periodically (de- 
loino) from the aircraft component. This requires the 
desion and installation of complex mechanical or thermal 
systems. Those systems are usually designed as an add-on 
or retrofit to an existing component. The second approach 
IS to design the component to eliminate or at least 
■intmize the adverse effects of ice accretion, such a 
component iould not alio, ice to Accrete, or the ice 
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methods are difficult to formulate, including all the 
necessary independent variables, and can not be used 
accurately to extrapolate beyond the available data base. 
An analytical approach does not suffer these limitations. 
Properly formulated, this method will not only reproduce 
existinq experimental data, but can be used to evaluate 
new airfoil deLigns. The theoretical model may also 
generate new insight into the icing problem. 

To be most useful an analytical jnethod must be as 
self containea as possible. That is, not rely on 
experimental results as input to the analysis. The 
analysis must contain in addition to the aerodynamic 
analysis, a model of the ice accretion process. The only 
inputs to the problem shou’d be the atmospheric iciJig 
envii'onment arid flight conditions of the airfoil. The 
logical first step in such ^ln analysis would be to 
initially study only rime ice accretion. Here the 
thermodynamics can bo ignored and the droplet dynamics and 
aerodynamics can be emphasized. Rime ice is streamlined 
in shape and conventional methods of aerodynamic analysis 
for- uiiseparated flows. ^.ce_a .p.p. 1 i .ca.b 1 e . Concentrating on . 
rime ice initially would provide insight into the problem 
while allowing time for the further development of methods 
-fjor dealino with icing thermodynamics an^the a nalysis o f 
separated flows. — 
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aircraft icing areas such as 

thermodynamics have rGcoivod recent attention, the 

analytical prediction of the aerodynamic performance^ of 

iced airfoils has not bwn .s tudied. Little experimental 

>fork has been done since 1953 and no attempt has been made 

to predict the performance degradation experienced by iced 

airfoils since Gray*s empirical method [3] of 1964. liie 

analytical prediction of the aerodynamic effects of ice 

accretion on airfoils is then an important gap in our 

knowledge of the icing problem, in a joint NAS^ and Yhh 

workshop on airrjLa.ft icing held at Lewis Research Center 

in 1973, the needs for new icino research were discussed. 

In his presentation Milton A. Beheim stated f4] 

... a renewed effort on icing effects on 
airfoils is needed not so much to 
refine ice protection systems as was done 
in the early 1950's but to determine the 
performance sensitivity to. ice accretion 
effects so that, airfoil selections can be 
made to avoid. using a protection system 
whenever possible. Particularly for 
general aviation applications it ..ay even 
be possible to evolve new airfoil 
geometries that minimize the 
possibilities of ice accretion and its 
deleterious effects on performance. 


This paper focuses upon the analytical treatment of 
two dimensional airfoils exposed -to a rime icing 
environment. New aircraft technology has generated 
requirements for an increased understanding of _the icing 
phenomena. This re -ex a tni nation of the icing "problem, this 


7 


tiinr with tho aid of hioh speed computers and modern 
numerical methods* promises the improvements in icing 
technology necessary to increase the utilization of 
general aviation aircraft and helicopters in adverse 
weather • 
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II. KKVIEW OF 

Tho qrowth ol ice on air 
df’croaso in performance and a 
the subject of scientific res 
tliif ortunately most of this wo 
year period which was conclud 
acio. Onlv now is there an at 
coordinate additional researc 
the proaress made bv early re 
need to continue this work, t 
literature is presented in a 

F,arly researchers disagr 
responsible for ice accretion 
the review of early work by B 
later French report [ 6 } in 19 
was performed in Germany and 
Countries. The u.S. made lim 
study- of aircraft iciflcj befor 
important hei-na the Oevelonme 
systems. These inflatable de 
and built by the P. F. Goodri 
on aircraft beoining in the 1 


LITERATURE 

craft components results in a 
5;aftey hazard w hich has been 
tjarch for over fifty years, 
rk was conducted in one ten 
eu almost twenty_-f ive years 
teii-pt to organize and 
h. In an attempt to clarify 
searchers, and document the 
his review of aircraft icing 
historical perspiKitive . 
ef'd on the physical phenomena 
. This is quite evident in 
leoker [5] in 1932 ana a 
3b. Much of this early work 
other Western Furopean 
itcd contributions to the 
e 1940; probably tho most 
nt of mechanical de-icing 
-icino boots were design€?d 
ch Company and wore installed 
930's. 


This same type of 
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systeui is still i-iv use touav. 

h major step forward was made in 1940 with the first 
mathematical fordiulation of water droplet trajectories, 
u, 1, Taylor [7] developed the differential equation 
-qovornir.q droplet trajectories £or the special cases of 
constant draq coefficient and Stokes law drag. 

Calculations were performed and the appropriate similarity 
uarametors extracted for a few simple cases. Taylor 
suggested a scheme for the numerical solution of the 
equation for more complicated cases ^such as the flow about 
an airfoil. This work was continued by Giauret [8] who 
performed the numerical solution of Taylor's equations by 
hand calculation. Giauret furthered the work of Taylor by 
combining droplet trajectories to determine the local mass 
flux on the airfoil surface and the total collection 
efficiencies. 

The publication of icing research in the open 
literature was discoirtinued daring the war years of 194 1 
to 1-945. However immediately after the war, .perhaps due to 
the need for all weather military aircraft made clear by 
the war, icing research flourished until the mid 1950 's. 

After the Second World War the United States* National 
Advisory Committee for Aeronautics, NACA, began an 
ambitious program. The bulk of this work was conducted at 


the NACA's lewis Rosearc* Center from 1945 to 1955. The 
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research was directed toward defining the natural icing 
environment, determining itK effects on rej'-rese ntat ive 
aircraft co-npujnents, and dosignim techniques for ice 
protection [4]. Great progress was made in understanding 
the icina procers and in protecting the aircraft from its 
hazards. The classic reference in the area of droplet 
trajectory calculations was published by Lanamuir and 
BTodgett [9j in 1946. Here the droplet trajectory 
equation is presented f or an arbitrary drag coefficient. 
The entire problem ol trajectory calculations is presents.! 
in a form similar to tJiat used today. Using a 
differential analyser the droplet trajectory about a 
cylin-de.r, sphere*, and riblKin were solved numerically and 
the collection efficiencies wex'e j-resentea for several 
cases. In addition, the modified inertia p>arameter was 
presented as a means to simplify the analysis by reducing 
the inertia p«jrameter and the free stream Reynolds number 
into a sinole dimensionless parameter. 

This method of numerically determinq droplet 
impinaement on aircraft components was used extensively by 
the NACA in the late 4C*s and early 50 *s. These 
rf’soarchors were greatly hampered- by the lack of high 
Pf-KH'd diaital computers and numerical solutions for the 
flow about an arbitrary bo<3y. As a result calculations 
were often made about bo.!ies for which the flowfield could 
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bp solvp>1 analytically. Droplet tra jectoric were 
calcnlatei? about cylinders [10,11], sptiores [12], and 
douKowsKi airfoils [13,14]. Arhitrary airfoil sections 
were first handled by Boiqrun [1b] usinq an empirical 
aperoa-cii based on droplet tra -jectories about JoukowsKi 
airfoils. nrun,aallaqlier ,and Vogt [16] used a vortex 
substitution method to generate the fl-owfield about an 
arbitrary airfoil. This approach required a wind tunnel 
tost to measure the surface velocities on the airfoil 
before the Vortex st_rennths could be determinoil . 

The method was used extensively by the NACA [17-19] 

I 

I to analyte the droplet, impinoement characteristics of 

i 

I 

j airfoils. Fxtensions of this? analysis wore mude by 

! Serafini [20] to a supersonic airfoil and by Dorsch and 

”run [21] to a swept winq. Droplet trajectory 
I calculations were also ierf.orJned about axisyniaielric bodies 

[22-Z4] to simulate the nose of an aircraft or aissle. 

I The trajectory calculations made by NACA researchers 

proved to be very accurate and provided valuable in: qht 

^ into aircraft icing, data for the design of qe-icing 

, system!', and, guidance to the experimentalists. 

Hally in the NACl. icinq proqram an extensive study was 
made of the natura] icing environment. Numerous 
experimental :;tudies were perforned to determine typical 

I 

com bin ati on s of cLcuix.. proper ties -s uch as horizontal and 

i 

L 
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vertical si7o, droplet diameter, liquiu water content, and 
air tempera tur«' experienced by aircraft. These data were 
compiled and rutiitnarized in three reports [25-27] which 
wore ultimately used to compile the FAl; Part 25 Appendix C 
[28] icinq envelope. This icino envelcfie is still in use 
and defines the ranac o^■ conditions over which any do- 
icing system must function to obtain FAA certification. 

Meuiy experimental studies were conducted in the NACA 
six by nine foot icing tunnel located at NASA Lewis. One 
important test prograO: developeu the dye-tracer technique 
for ex ner i men tally obtaiHing impungement characteristics 
of arnitary bodies [29]. In the dye-tracer technique a 
body is confiqured with blotter paoer and exposed. to an 
airstreain containina a dyod-water spray cloud, Tne 
blotter paper is then calorimetrically examined in order 
to obtain local collection efficiencies, total collection 
efficiencies, and maximum rearward extents of impingement. 
This technigut! has boeii used on airfoils [30^31] and other 
bodies [ 32,33 ] and provides the only direct experimental 
data for use in the verification of droplet trajectory 
calcul ation s. 

Airfoil icing oxj'eriments conducted in the icinq wind 
tunnel served two main objectives. These tests docuuti-cnted 
the Change in airfoil performance char acteris tics due to 
ice accretion while also serving as test beds for new de- 
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icinq and anti-icinq sy:rtem?. In the first test [34,35 ] 
no quantitative measure was made of the ice growth. 
Aercul^jLoaJiijCL data was obtained from a heated wake survey 
probe measuring the changes in drag, while lift and moment 
—coefficient changes were not measured. The tests were 
primarily to evaluate 1 ice protection systems. Bowden 
[36] in 1456 presented a fairly complete aerodynamic 
evaluation of icing effects on a NACA 0011 airfoil. A six 
component force balance system vas- used to enable the 
measurement of changes in lift, drag, and pitching, moment . 
As in earlier tests the documentation of the ice shapes 
from which the aerodynamic penalties resulted was only 
described aualitatively . 

The most complete airfoil icing analysis performed is 
reported by Gray [1,37]. Here theoretical and 
experimental i T.pingeB.on t efficiencies, ice shape 
measurements, and an aerodynamic analysis was performed on 
an -NACA 65A004 airfoil section. The experimental and 
theoretical impingement characteristics compared well for 
some cases,, but the failure of the predicted values in 
some situations, was not understood. Gray [34] presented 
the first empirical relation to be used to predict the 
chanaes in drag coefficient due to icinq. This equation 
was based on the NACA 65A004 icing data and was good only 
for this particular airfoil. 
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In approx iiTiatcly 19b0 icino research at the NACA xas 
stopped. The advent ot jet ennined aircraft reduced the 
icino hazard and required that research efforts be shifted 
to new areas. At its coit'pletion the NACA prourarr had 
provided oocd ice protection for the aircraft of the day- 
The analytical prediction of iirpinqement rates had begun, 
but no inothods for ice shape calculation or the resul ting 
airfoil performance degradation were developed. 

Experimental results were confined to only a few 
soecialized airfoils and had consisted primarily of ice 
protection sy5;te;n evaluation. Two compilations of NACA 
data were published in 1964 , Gray [3] compilec^ all the 
iced airfoil drag data to expand his empirical equation 
and i'owden et . al. [ 3B ] present^id an exhaustive summary of 
existing aircraft icing technoltxj y , 

Interest in aircraft icing research was renewed in the 
early seventies in Europe and Canada. These studies have 
primarily involves with the thermodynamics of* the ice 
accretion process. LozowsKi, Stallabrasr, and hearty [39] 
in 1d7u presented a summary of t-hermolynamic modeliing and 
their current state— of ~the-art approacli. All of these 
studies are hampered by the lacK of good droplet 
impinoement methods. Research has been conducted in 
Western Europe in several areas which are summarized in 
reference 40. Recent aerodynamic studies have been 
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conducted in Sweden and the Soviet Union [ 2 ] to determine 
ejcrerimentally the performance of iced airfoils. Similar 
tests conduct Cvl in the United States have in general not 
been published since they wore conducted by manufacturers 
for icing certification and not government sponsored. One 
recent exception is the work by Wilder [41] from Boeing, 
wilder presents the results of theoretical impingement 
calculations, experimental ice shape correlations, and 
iced airfoil tests. Unfortunately little information is 
provided as to the analytical or experimental methods used 
to obtain these data. 

Recognising the need for an organized, icing research 
effort in the United States, NfiSA Lewis Research Center- 
established a program of icing research in 1980. The nasa 
program includes a broad range of research objectives. 

The evaluation of de-icing systems and anti-iedng systems 
[42] has recently begun in the Lewis Icing Tunnel. 
Analytical efforts include a three dimensional droplet 
trajectory code [43], and preliminary results of this 
-d-issertation [44]. Hopefully the need to apply current 
technology to the icing problem, as revealeg-by this. 
review of past research efforts, will be met by the 
current N-frS-A icing research program. 
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III. THEORETICAL ANALYSIS 


Before th(^ aero^lyne mic t>eif ormance of an airfoil with 
rime ice can he determined, the geometry of the ice 
accretion must first bo calculated. This section presents 
the theoreti cal method for the prediction of rime ice 
shapes which accrete on unprotected airfoils. Therefore 
tne first ste] in the theoretical analysis is to formulate 
and analY 2 C the equation governing the trajectory of a 

sinule spherical particle in a moving, fluid,..., , , 

Trajectory Equation 

Aircraft rime icing occurs when super cooled water 
droplets im;>act the leading edg«i of an aircraft component. 
There- droplet s have cia^ieters of 10 to hU microns [20] and 
pxoerience Reynolds numlxprs low enough to ensure that the 
particles reniain spherical in shape [43]. For rime icing 
clouds the liquid water content which exists rarely 
exceedes 2.0 (jrams of water per cubic meter of air. Due 
to this low concentration of water droplets in the free 
stream, the flow may be considered uncoupled [4lj]-and the 
influence of the droplet on th€» flowfield ignored. 

By applyina Newton's Secemd Law, P-ma, to the 
particle, the governing equation can be derived. This 
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equation has been piesenterl by Soo [46] anil Rudinger [47] 
as “ 



This equation way be significantly reduced for the present 
application. For a water droplet moving in air the 
density of the particle is much greater than that of the 
fluid. Therefore, the pressure gradient term, P, and, Ma, 
the apparent mass term may be neglected [46,47]. The 
foUE-Ui. term in equation (1) , B, represents the Bassett 
force. This term accounts for .the devi ation of the flow 
pattern around the particle from that of steady state and 
represents the effect of the history of the motion on the 
instantaneous force [47], This term is significant if the 
particle density is of the same order as that of the 
fluid, or if the particle experiences large accelerations- 
Droplets can experience large accelerations when in the 
leading edge region of an airfoil- Morment [43], using 
™the work of KeinL_[ 48 ] and Crow [49], has shown that for 
the icing problem the accelerations experienced by the 
droplets are not large enough for the Bassett term to be 
significant. Therefore the Bassett term, B, can also be 
dropped from the analysis. With these assumptions, eg. 

(1) reduces to. 


2 _ 

°‘(^) = B + mg 


( 2 ) 
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For the_ S!Tiail water dropletr. considered in the icinq 
problem, and for the small time scales involved, the 
aravity or set^tlino term may in 9 eneral be dropped from 
the analysis* Fiowever, it will be retained here to allow 
a more qeneral application of the method* The viscous 
drao term, D, can be expressed in the conventional manner 
as 


B ■= i 





Where S is the cross sectional area of the sphere and Cq 
the d.Tag coefficient derived from experimental results. 
Note that here the draa is evaluated using the slip 
velocity, that is the velocity between the droplet and the 
lcx:al airstream. Substituting in for the drag and 
dividing through by the mass, eg. (2) becomes 


. 3 !l£a|a.dx|/- g 
dt^- 4 6o I dtl\ dt/ 


Nondimensionalizirig eg. (3) yields 


n ■ K h I 


(3) 


(4) 


which is the governing equation for a droplet trajectory. 
The nondinrensionalizaticii was performed with respect to 
the characteristic velocity, U, the free stream velocity, 
and the chTaracteristic length, c, the oirfoil chord. 

Dif ferontiation is with respect to nondimonsional time,! , 


19 


whf>re t = Ut/c, 

In eq . (4) tifo d iaensionless parameters occur. The 

inertia parameter, K , is 


K = 

18cp 


(5) 


and is essentially a nondimensional particle mass. The 
second parameter# Pr. is the Proude number 




, U 

Vcg* 


( 6 ) 


which is the ratio of inertia to gravitational forces. A 
third similarity parameter appears due to the form of the 
CdR/24 term in eg. (h) . 

The drag coefficient of a sphere in a oon^ccelerating 
stream has been measured as a function of Reynolds number 
by many researchers. Sphere drag is also in general a 
function of particle BAriti. number. However, for rime icing 
which occurs at low flight velocities, the particle Kach 
numbers are low and the compressibility effects on sphere- 
drag are not significant. Here Reynolds number is based 
on droplet diameter and the relative velocity between the 
stream and the particle. Reynolds number as used here is 
given by 


= P*sulu-nl 
y 




A standard drag curve has been established from these 
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results and is presented in Schli.Qhti.nq [50]. Por low 
Reynolds numbers the well known classical Stokes solution 
for sphere dran is 



However this thearetical result is for creeping motion and 
is not valid for the hiqher Reynolds numbers experienced 
by icinq droplets, stokes drag law is however a limiting 
case used to establish empirical fits to the standard 
sphere drao curve good for hiqher Reynolds numbers. 
Langmuir [9] presented cme of the earliest empirical fits 
of the sta.ndard sphere drag curve given by 

-^ = 1 + 0,197R°'^^ + 2.6 X 10'^R^-38 <8) 

This equation provides qood drag coefficients up to a 
Reynolds number of 1000. A somewhat simpler form proposed 
independently by Klyacho [51] in 1934 and Putnam [52] in 
1961 is 

1 + i r2/3 <9, 

24 6 

Both eq . (8) and (9) represent good fits to the 

experimental results as do several other similar equations 
proposed by other researchers. The standard drag curve, 
StoKes law, eq. (8) , eq . (9), .and the recent and— more 
accurate measurements of Beard and Pruppacher [53] are 
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compared in Table 2. 

Table 2 Comparison of Particle Drag Coefficients, Cp 


R 

Std 

Stokes Eg. (8) 

Eg. (9) 

Beard and 
Pruppacher 

0.01 

2400. 

2400, 2426. 

2420. 


0.1 

24 3. 

240. 251. 

249. 

242.7 

1. 

26 .9 

24.0 24.5 

28.0 

26.45 

10, 

4.33 

2.40 4.43 

4.26 

4.149 

100. 

1.09 

0.24 1.14 

1.10 

1.073 

500. 

0,568 

0.048 0.588 

0.552 


1000, 

0.469 

0.024 0.477 

0.4 24 


The empirical 

fits. for this sphere drag coefficient 

including 

equations (8) and (9) are 

of the general form 



N Yi 

— = E C.rYi 


(10) 



24 i=l ^ 



Osing eg. 

(7) for 

P , eg . (10) can be written 

t as 



N ■ _ . 
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|Yi 

(11) 

where Ry 

is the free stream droplet 

Reynolds number 



1, P<SU 





V — 


(12) 

Therefore 

f since 

the droplet drag coefficient can be 

expressed 

in the 

form of eg . (1 1) , 

the Stokes parameter. 

CdR/24, appearing 

in eg. (4), yiejLds the _addi ti onal 

similarity parameter Ry, the free stream droplet Reynolds 

number . 





The trajectory of a liquid droplet, for 

the rime icing 


problem, has been shovn to be governed by the differential 
equation (4). Eg. (4) contains the three similarity 


parameters Ry, K, and Pr. In the next section Ry and K 
a Ee„. combined into a single parameter which greatly 
simplifies the analysis. The flowfield which enters 
equation -(4) as u, will also be discussed in a later 
section. 

Trajectory Similarity Analysis 

In the derivation and discussion of eg. (4) it has 
been shown that the droplet trajectory, ignoring the 
flowfield, is a function only of the three similarity - 
par-ameters Ry, K, and Fr, and the initial droplet 
conditions. To simplify this an. -^ysis the Proude number, 
Fr, will be dropped, since it can be shown to be 
negligible for the rime icing problem. The scaling of the 
gravity force and other terms in eg. (1) will be discussed 
later. Nealcctinq the <iravity term eg. (4) becomes 



Now the, tra jectory depends only on F and K, assuming the 
initial conditions in naidimetisional form are con.sta.nt. 

The identification of the proper similarity parameters 
for a problem is very important. Not only do the 
parameters simplify the analysis, but they also aid in the 
presentation of expe^rim-ntal and numerical data, and serve 
as -scaling parameters in the design of scale model tests. 
Por-aircraft icing scale model- tests, using Rtj and K to 
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establish test conditions violates other similarity 
parameters. For example, often only the model speed and 
droplet diameter can be varied. ?!olding and K constant 
then requires that 

i„-(Sa)s and 

hs a result, for small scale models the test velociti.es_ 
are very large and violate the Bach number scaling of the 
flovjfield. Similar problems in the scaling. of drops in 
aircraft wakes have been reported by Ormsbee and Bragg 
[54]. Recent icing tests by a Swedish-Soviet research 
group [2] chose to ignore the Reynolds number scaling and 
hold only the inertia parameter constant in an attempt to 
avoid this problem. 

Rethods are available to alleviate this scaling 
problem by reducing the number of similarity paraneters. 
Coiiujininq the similarity parameters Ry and K into one 
parameter would also greatly simplify data presentation. 
The first attempt to combine and K was made by Langmuir 
when he presented the modified inertia parameter, Kq • 

This parameter will be discussed here and a new derivation 
presented which for the first time yields an analytical 
solution. In addition, a method is presented which is 
much simpler to use, and in many cases mere accurate. 


Inert i a PaiciBu^tor : Tht> w odi iied inertia 

ii!TU’t«^r , Kq» was prt"Sonti>d by Lanqmuir [9] in 194b to be 
us(m 3 to-i>resen1 aircraft icinq.data. In fact, this 
paramrti'r is still in wide use in the aircraft icinq 
com'nunity athoiwh no thf\ic.otical proof of its validity is 
available [ 3f’ J. Currontly no clos*-*d form solution for the 
parainetpr exists and a qiaphical technique or curve fit to 
the numerically uenerated data is used. Here the 
parameter will derived from the qoverninq differential 
equation and a dosed form solution obtained. 

The modified inertia parameter, K^, is defined as 


K = K 
o I 1 , 


where K is the inertia parameter and \/\g is the ratio of 
the traqeciory of a tiioplot in still air, with an initial 
Reynolds number of Ry and qravity neqlected, divided by 
t-he same trajectory of the droplet if the drag is assumed 
to obey Stokes law. So Kq combines K and l<y into a single 
parameter since \/x is a function of Ry alone. Lanqauir 


showed that \/ is given bv 


^ A. r'^ 

J 


Using the statid.ird sphere draq curve lor Cdii/24 , Lanqmuir 


performed this intearation numerically t6>- qenerate X/Xg as 
a function of Rn which is still in use today. 
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By using the differential eq . (13) we can exanine 
more carefully the origin of Langmuir's Kq parameter. It 
is not clear from reference 9 if Langmuir derived Kq in 
this way, but tho_lxasic relationship between Kq and the 
governing differential equation was suggested in 1952 in 
reference 55. 


By rearranging eq . (13), it becomes 

q =(G-n) 

. 24 / 


(16) 


Here CdR/24 is a complex function cM. R, with R varying 
along the particle trajectory,. If some suitable -average 
of the term on the left hand side- of eq. (16) could be 
found, the Ry and K parameters can be combined into a 
single similarity parameter. Assume that the particle 
experiences Reynolds numbers from zero to R the value 
based on the free stream velocity. Then averaging this 
term yields 


K = 

R. 


J % 


U dR 
24 


(17) 


The modified inertia parameter is merely the average value 
of the single coefficient which appears in the droplet 
trajectory equation (16). is not an exact similarity 

paraJueter, but does have valid the o-retica t justification 
as it is a straightforward simplification of the governing 
particle trajectory equation. The modified inertia 
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I'draiiu'ter proviuot. good data correlation, provided the 
rang*’ of Reynolds numbers experienced is consistent with 
the range zero to R^. 

A closed t orm solution for Kq can be found if an 
integrable torm of the droplet drag coefficient is used in 
oq. (8) . Putiia? and Klyacho [51] devoloiied such an 

equation valid ui ".’nol ds number of 1000 as 

‘‘'7i + 1r2/3 (18) 

28 6 . . 


Following the work ol Putnam, and after considerable 
•integration and algebra, a closed fora of Kq is given as 



This equation is within one percent of Langmuir's 
calculated values until Ry approaches 1000 where 
Lanqmuir's values deviate from those of eg. (ig). This is 
due in part to the different droplet drag values used, 
eg. (3) and (Id) , and probably same accumulation of error 
in the nume-rical procedure. 

The lower limit of nj . (1*)) can be used' to check the 
derivation of K^. Py definition must appr^aach the 
inertia parameter for small values of the PeynoLds nuiber 
where the particle drag is essentially governed by Stokes 
law. By expanding the inverse tangent function and taking 
the limit as R„ approaches zero, eg. (19) reduces as 
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to Kq equals K. By examining eg. (19) as Ry 
approachep infinity Kq takes the form 

• =" Kq = 18K <20) 

It is also interestinu to compare the curve fit develo ped 
by Stallabrass and Lozovski [39] for Kq where 

K = ZTSTT-, (21) 

° 1+0.0967 

This compares well to eg. (19) ; note the similarity in the 

6367 -2/3 

1/Ry' term in eg (21) and the expression in eg. 

(19) . 

The use of eg. (19) should improve the usefullness of 
the existing icina data correlated with Kq . Eliminating 
interpolation or curve fits to Langmuir's tabulated data 
should also improve accuracy, Eq . (19) could be used to 
reduce other droplet trajectory data, however, the 
analysis to follow will result in a parameter which is 
easier to use and- more accurate and versatile than the Kq 

parameter. 

Trajectory Scaling Parameter : An alternative 

approach can be taken to simplify the single coefficent 
appearing in the trajecto^ eg. (16). Instead of assuming 
that CdR/24 is a constant, as wu.s done to derive Kq , here 
assume that 


CdR 

24 


CR^ 


( 22 ) 


which is tho first t^^rm of the general equation (1 ). 

This appeals us a struiqht lino on the log-log plot ol 
CdR/2« vs. R, liouto 3. A-similur approxi mation _has been 
made before by nrms bc'p-aiul Bragg [54,56] and by Ar«and et, 
al. [57] to scale droplet trajectories. The txajectoiy 
equation becomes 


/CR 


U 


= 


V K 


^ '^ju-nr(u-n)- 


(23) 


Now define the trajectory similarity parameter, K, as 

(24) 


K = JL. 


where the coe^fficient in eg. (23) has been inverted to 
follow the convention established by the mcKiified inertia 


para iiiC ter • 

The appearance of the |u -ppterm in eq . (21) 

simplifies the use of this parameter while decreasing the 
expected incTrease in accuracy over the parameter. 

Since a y occurs outside of the K term, C and y cannot in 
qeneral be functions of but must be chosen from a 

sinule best fit of CdF/24 - CK^ over the entire range of 
Reynolds numbers to be experienced by all particles under 
consideration. Then after C and y have ber^n chosen for a 
particular application, a simple parameter combining K and 
R is available to be useif for data presentation or 
establir.hing scale model test conditions. Note that it 
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thf oravity term need be included, this re<]uires only that 
the Proude jiumber, Fr, also be considered in addition to 

K . 

A careful analysis of the nodified inertia paraaeter, 
Kq, and the trajectory scaling parameter, K, shows that 
th«? two parameters are related. Lf the appro xiaate drag 
law of eq . (2 2) is assumed and used in eq . (17) the result 
is 

K = ^ r 

c(i-y)r^7 

For this special case then differs from R by only a 
constant. While the general form of given by eq . (19) 
is more complicated, it too can be seen to be in a 
functional form similar to that of K. Taking the limit of 
both Kp and K as Ry approaches zero yields just K in both 
instances. As Ry approaches infinity the limit of is 

Kq = 18KR-2/3 

as given earlier in eg. (20). This is exactly eq . (22) 
for K if ^ = -2/3 and 1/r = 18. Therefore it has been 
shown that and K ere the same within a constant if a 
simple drao law is assumed in deriving Kq. So Kq and K 
are certainly closely related but do vary slightly in 
their workino ranuc. 
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A irct hoil for r><:*loctiua the Y to be used in the K 
par a me tor inurt ix' dc'lt»mihed • The selection of a Y 
dopcndn on the ranoo of Peynoldi; numbers to be experienced 
b^y t-he uarticlr^; durincj their trajectories. Since a Y 
term occurs i li the datferential equation outside the K 
t('rm (see 23) unlv one y may be selected for each 

arn^licaticTn of h. So for a scaiino applicatic»ri a 
different Y aiav bo selected for each particle considered, 
however, when K is uied. to present data, an avoraqo valiu' 
of ') must be used which is qood over all uossible particle 
traiectories to be presented. The value of the constant C 
appoarinq in the K pardiiiotor has no effect on the use of 
the term and is considervxl to be equal to one for 
con veniencr^ . 

St^foro Y ran be determined the 'loyrrolds number ranqe 
experienced by a rime icinq nartrcle must be estimated. 

The suner ctKjled water droplets are assumed to be at rest 
with respect to the atiDosphoro initially, and therefore 
the lower bound on the Ueynolds mrinber ranoo is zero. The 
d-ro^Htet^exper iences its maximum Hcvnolds number when it is 
in the immediate proximity of the airfoil leadina eduo 
where the voiocitv (iradients, and therefore th%*' relative 
Velocity hetwt'f'n the dr<.»i anil the fluid, are laree. In 
fiqurr 4 the Heynoldr. nuui!)ors experienced bv droplets 
alona their traiectory is shown. This i iiLor nia tion was 
qent*ratP‘1 urinu the ('emeuter code described* in Section IV. 
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Th*^ ilroplpts wore starttvi five c ho rds in front of a 
typical qoneral aviation airfoil opcratiiiq at a cruise 
condition. Note that all the particles experience 
Reynolds numbers in the StoKes law range for the first 
nint ey percent ol their trajectories. Only as the 
droplets approach the btx1y do the Reynolds numbers 
increase dramatically. This analysis has shown that the 
droplets usually experience maximum Reynolds numbers of 
less than one -ha If Pjj . 

Usinq this; information on the typical Reynolds number 
range alonq with fiaure b a value of y can be determined. 
Fiqure 5 summari7es the results of a least squares fit 
program which calculates the value ol Y which provides the 
best fit of the approximate sphere drag expression of eq . 
(22) to the standard sphere drag curve. The fit is 
performed from a Re>ynolds number of xero to R . It has 
been found theit for the aircraft icing problem a Y of 0.35 
represents a good average value to be used for preliminary 
scaling calculations and for data presentation. To select 
a "I to use in scalinu a particular droplet, the average 
value of R for the full scale and scaled particle is 
found and then figure h tran be used to deterinine y. Xn 
general this is an iterative procedure, hut by using 
1 = 0.35 to select the initial scaled hy it converoes 
rapid ly; ur:ually thf^ lirst st*^p is sniticivntly accurate. 
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_ » i^ystt‘nio t ic procxniuro has presented to reduce by 

ono thtr‘ numlK^r of sindlarity paraniotors yovprning this 
class of part icle trajectories. The ui<'tho(j of La-nqniuir, 
previously^ little understood, has been derived from tne 
qoverninci dillerential eHination and a clost*a form solution 
has been nreronted. This r^s^Ailt should clarify the 
theoreti-Ci^-l — basis for the modified iiu.'rtia parameter anu 
iriaKe the existiiKj data correlated usinq — easier to 
i nt er^ ret . 

^ new dimensionless number, K, the trajectory scaliria- 
parameter is derived. This pciranu^tor is more accurate and 
versatile than the modified inertia parametei. The 
trajectory scaiina parameter may bo used to simplify any 
trajectorv analysis. All that is roauired is the 
dot ern.inc tion of the oxwnent, >| , in the approximate drao 
law u^ed in derivind K. The exponent u.av hk found by the 
f ol low inn procedure : 

1} Determine the ranue of Keynolus numbers 
oxiH^rienced by the class of particles for 
which the K {'arametor is to be used. 

2) T^v us*ino a least s<}uares or other best fit 
rcheme, determine the y for which the 
approximate drap law hert. fits the standa.rd 
drao curve in th*.‘ Reynolds number ranae 
of interest. 

^xnerinien tal and numerio^i re.^ults in support of th<^ K 
farameter, anrt a comparison of and K will l>e proserrte^ 


in Sort ion VI 
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Flowfigld ■ — 

To calculate the trajectory of a particle in the 
vicinlty—ol an airfoil the detailed flowfield must first 
bo determined . The dimensionless flowfield velocity 
appears in the differential equation (4) as "u, and also in 
the Reynolds number, P. The effect of a compressible 
flowfield on water droplet trajectories has been studied 
[10] and found to be negligible for cases up to the 
critical Hach number. In addition, the viscous effects 
near the leading edge of an airfoil are confined to a very 
thin boundary layer. Since for most applications the 
water dropletr* only impact the airfoil near t!ie leading 
eduG, the effects of the viscous region near the airfoil 
are assumed negligible. It is therefore sufficient for 
this purpose to describe the flowfield about the airfoil 
by an inviscid, incompressible, potential flow solution. 

Both singularity and conformal mapping methods are 
currently in use for predicting the flowfield about an 
airfoil. Both methods were used in some form by the NACA 
to maKe droplet calculations i.n. the 1960*s. This present 
analysis uses a mouified version of a transformation 
scheme for arbitrary airfoils first presented by 
Tneod-errsen thP,f>9]. This method as formulated by Woan 
[ftO], replaces the JouKowski transformation used by 
Theodorsen for the first step by a Karaan-Treff tz 
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ti'anr:f orniat ion . This provides a better near circle for 
airfoils with finite t raili n'? edgji_a.n.cLifi.s . Tlie second 
important feature of this motnod is in the solution for 
the exact circle. For most applications only tne velocity 
on th»‘ airfoil surface is desired. To. oi»i:.ain only_±be 
surface velocities a simplified approach is available 
wliich eliminates the need to calculate all the 
coefficients in tbe—com plex Fourier series which 
transforms- the near circle to an exact circle. However, 
the Fouri(*T coefficients are needed to dete.rniino- the 
velocity in the flowfield at some point not on the 
airfoil. The analysis of Woan provides for the 
determination of sufficient Fourier coefficients to solve 
for the entire flowfiel>3. 

The TheodcjirtMi methal has been criticized in the past 
for tailinu to reproduce the flowfield accurately, 
particularly in the vicinity of the leadinc edqe. It now 
appears that this was a characteristic of early numerical 
schemes and net of thf* method itself. The velocity 
distri butioii computed about a typical general aviation 
airfoil by this method is compared to the sophisticated 
sinaularity method of I ppler [61] in fipure 6. Only the 
leauinu edoe re.jion is stwjwn where the calculated velocity 
distributions are practically idejvtical. The comparison 
of the surface velocities over the rest, of the airfoil is 
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also excellent. This denionstrates the validity of the 
Theodorsen flowfield code. 

The method provides an accurate velocity anywhere in 
the flowfield . This information i s use d in the solution 
of eq . (4) . In addition this method has proven to be very 
successful in handling leading edge shapes required later 
in the analvsis. 

Droplet Impingement Parameter . 

By analyzing the information gathered from several 
droplet trajectory calculations much useful information 
can be extrapolated. Glauert [B] first combined droplet 
trajectories to determine the mass of water strikina a 
circular cylinder as a function of theta, the angle 
measured from the staonation point. Langmuir [9] extended 
Glauert's analysis to determine the inass striking an 
arbitrary body as a function of S, the arc length along 
the surface. The analysis presented here will follow that 
of Lanamuir with some extensions, particularly in the area 
of clouds containing distribution s_of droplet K iyes. 

kssuming the droplet trajectory information is 
available, the first step is the calculation of 0, the 
impingement efficiency. The impingement efficiency is a 
dimensionless mass flux of the material impinging at a 
particular point on the airfoil surface. B is 
nondjjiie,as.i.oaa.lized wi th respect to the mass flux in the 
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frt'H stream. An impingement efficiency of one is just 
that in the free stream, o r it is the dimensionless mass 
flux on an imaginary flat plate, which does not alter the 
free stream flow, placed perpendicular to the free stream. 

The impingement efficiency on an airfoil surface can 
be deduced from figure 7. The position on the-airfoil 
sur-ftrce is given by S, the arc length along the surface 
measured from the leading edge. S is measured in chords 
and is positive on the upper surface and negative on the 
lower surface. The vertical position, dimensionless with 
c, in a plane perpendicular to the free stream is given by 
y o* The mass of water droplets between the two particle 
trajectories a distance apart in the free stream is 
distributed over a length 6s on the airfoil surface. As 
the length 65 approaches zero, the local impingement 
efficiency becomes 


dyp 

dS 


(25) 


Note that in the free stream dy^ equals 6S so that 
B - 1 as .equired. g can now be calculated by taking the 
derivative of the as a function of S curve derived from 
individual droplet trajectory calculations. 

The total mass flow rate of water caught per unit 
span by the airfoil is then given by 

M = UXc ^ ^ gds 

Sl 


( 26 ) 
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Here the limits S .. and S are, respectively, the maxiaua 

U Jj 

liisits of droplet impingenent on the airfoil upper and 
lower surfaces. 9y substituting eg. (25) for 0 in eq . 
(26) , n becomes 


M = UXcAyf, 


(27) 


The total mass collected by the airfoil then depends on 
A Yq, the distance in the free stream between the upper 
and lower tangent trajectories, figure 7. It is 
convenient to define an overall collection efficiency, E, 
to evaluate and compare the impingement or catch rates XLf 
various airfoils. The collection efficiency is defined by 
the rate of mass caught divided by that of the free stream 



Here h cart have two different values. Some researchers 
take h as the maximum airfoil thickness to chord ratio, 
while others use the maxinium projected frontal height, 
which is a function of angle of attack. This paper uses 
the later definition unless otherwise specified. 

The preceding discussion describes the calculation of 8 
where the icing cloud contains only a single droplet si2e. 
In general clouds contain a distribution of particle sizes 
about some volume mean diameter, VRD. To represent the 
total impinqement efficiency, B^,for a point ori the 


airfoil including the particle size distribution effect, 


the equation is 


/ ^max . 

\cl6/ 


(29) 


•^min 


Here 3((S,5) is the impinoemen t efficiency at a point S on 
the airfoil surface due to a particle size 6. Langnuiir 
[ 9]_has -defined four particle size distributions about the 
VMP whiah..are fairly representative of actual icing 
clouds. The distrib-u-tiuis are_gefined by V, the 
cummulati ve voluo’e of water in the cloud, as a function of 6, 
the droplet diameter. The (dv/d6) term in eg, (29) is the 
derivative of this curve and is a function of only 6. 
Considerinn the entire ranue of droplet sizer also 
complicates the calculation ol the total mass and 
collection efficiency, the total mass becomes 


M, 


.5. 


= UXc/ 6(d,S)(^)d<5dS 

'min ^ 

and the collection efficiency is 


(30) 


3lax 


_''y.o(«)(f V* 


(31) 


‘min 


Vdd/ 


Here Ay^Cd) in the difference in the y^ values in the free 
stream between the tangent tfagectories for a particle of 
size (S . The value Ay^ (6) can be determined directly from 


the analysis or is given by 

S 


Ayo (<5)= / 3 ( 5 . S) dS 


(32) 


U 
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The im piny ament efficiency, 6, and as a result the 
total mass caught, M, and the collection efficiency, E, 
can now be deterininecl by combining the results of several 
droplet trajectory calc ula tions. For the rime ice case, 
knowing 3 as a function of S and the free stream 
conditions permits the prediction of an ice shape. 

Ice Shape Calculation 

Using the information £Tovided by the 6 curve an ice 
shape can bo predicted for the case of dry accretion (rime 
ice). Glauret [8] recognized this relationship between 
and the rime ice deposit. However he was only able to 
give a pictorial representation of the shape by measuring 
out from the surface a distance prjQportional (to an 
arbitrary scale) to the local rate of droplet isipingement. 
wilder [^1] has calculated rime ice shapes assuming the 
ice grows out -normal to the airfoil surface, but has 
ignored the local curvatuire of the airfoil surface. Here 
the equation for ice growth will be derived including the 
effect of surface curvature and an arbitrary direction of 
ioe growth. 

Consider an area dA perpendicular to the free stream 
velocity vector. The mass of water passing through this 
area in a time At is 


m = U\At3dA 


^0 


Not. 3 th. collection .efficiency on the surface 3A . 
T|,p .olilJU! Ot ice-, >1A', rerr«s.-nt»J b, » is 


uxstMA 

Pice 


Koarranqiny aa>1 naaai.«.^ionaU^in=! 


I = A^p 


(33) 


.here Ac is a new similarity parameter given as 


A = 


(34) 


'ice 


r- nar he in ternrettc 3 as the 

The accumulation parauiet, 

in airfoil chords that would 
length of the ice growth in airi 

occar on an ioaoinary flat ,.late place porpen.Ucnlar to 
the free ntrea. in . tiu.. At. bote that 6= 1 on this 

T +Tnn naran^^t€r Qi^vorni^ ih^ riwo 
flat plate. The accumulation parameter m 

, Q r^it-ve has been determined. It is 
icing process once a g curve na.. 

the cross sectional area of an ice 
convenient to represent the croo. 

^hape in terms of Ac using the expression 


A - / '• 


AjngdS 




p..rton«ing Iho inloration the area boco.os 


A = A^AYo 


(3b) 


Since AC am Ayo -e both liacnsicnlenn , the atea P t.en by 
eg. (3b) has units of square chords. 
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Now usinq the concepts of accumalation parameter, Ac, 

and local impinoement efficiency, g, the ice shape can be 

determined. Figure 8 shows the ice growth (cross- hatched) 

on a small segment of the curved airfoil surface ds. Here (j> 

is the assumed direction of ice growth and r' the 

effective radius of curvature of the surface. (The 

effective radius of curvature will be defined later.) 

From geometry and noting that the ice area must equal he gdS, 

«2 

A + ^ = AgL. (36) 

2r' cos4» 

This may be solved for i, and is the general expression 
for the length of the time ice accretion, for a given Ac, 
at a point, S, on the airfoil surface. Here g, ip, eu)d r* 
are all functions of S. Two special cases of eg. (36) are 
of particular importance. 

The first case is to allow the ice to grow out normal 
to the surface. Here 4> = 0 and r' is just r, the loccil 
radius of curvature of the airfoil surface. Eg. (36) then 
becomes for- normal growth 

^2 

^ = AgS (37) 

Here a nonlinear term arises due to the radius of 
curvature of the airfoil, r. This term has been dropped 
by other researchers when calculating £. This assumption 
is justifiable for small values of Ac or for airfoils with 
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Note that when r iS lar«ie» 

a larc,e leadii>g edqo ra xus. 

V <= that the length of tne 
(37, sho-s^ tha can'easily 

^ nf the 

The i.portance area of the ice 

evalualei hy coaparing 

shape to the heeh sag 9 estea la 

‘ hath out alohg the particle tralector, 

.hich the ice ptoas ^ ,o the 

- (373. 1« ,36,. ,-e P iS 
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Two modps of ice growth, normal and tangent, have been 
discussed in relation to the solution eg. (36). However 
eg. (36) can be used for any ice growth s cheme i f the 
trajectory tangent in figure 9 is replaced by the assumed' 
direction of growth and r' and <J) are determined 
accordingly. No matter what scheme is used, after eg. 

(36) is solved for £ , it is easy to calculate the ice 
-s-hape by moving ovit from the airfoil surface a distance 
in the f direction. 

Time- -M fects 

As the ice accumulation builds on the leading edgje. of 
arr airfoil, tlie tlowfield must slowly adjust to the new 
boundary conditions imposed by the change in shape. This 
chanqe in the airfoil shape, and the resultinq change in 
the flowfield, will naturally alter the impingement rates 
on the surface. As tho impinaeiront rates chanqe, the 
shape of the resulting ice accretion will also change with 
time. Therefore the ice accretion process is a function 
of time, and, must be modelled accordingly if accurate 
analytical predictions are to be realized. The failure of 
initial icing rate calculations [37], or shapes based on 
them, to accurately predict the experimental results 
reinfoTces the need to include time dependence in the 
model. One method of modellino the effect of time is a 
time-=ste.pping approach. The time-stepping method assumes 
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that thp ico accretion can be broken down into a seri es of 
steady state processe.s. The accuracy of the method is due 
in part to the step size chosen. 

The scheme used to perform the time steppinq is itself 
relatively straiyht forward. bach time step csin be broken 
down into thfoe .par.ts: 

1) The flowfielu is generated 

2) The g curve is calculated from the pia.rticle 
tra jectoricK 

3) An ico shape is generated 

These steps are then repeated until the desired icing time 
is reached. li; 'ractico the procedure may be very 
difficult since the iced airfoil coordinates qenerated in 
step 3 may be too "rough” to permit the calculation of a 
flowfield. A scheme for smoothinn these coordinates is 
discussed in Section TV along with the numerical 
formulation of the pioblom. 
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IV. NUKERICAL PORWULATION 

The theoretical analysis presented in ii.ection 111 has 
been proqra^imed for commuter solution. This section 
describes the numerical procedures and computer codes used 
to predict tne rime ice growth on airfoils. The solution 
is formulated into three steps which utilize four computer 
programs. The three steps are: 

1) Droplet trajectory calculation including 
flowfield generation and the determination of 
impingement rates 

2) Kime ice shape calculation 

3) Iteration and coordinate smoothing 

Step 1 contains two computer programs, while step 2 and 3 
contain one each. A flowchart for the entire rime ice 
methodology is given in figure 10. Only the flowfield 
code was not written especially for this study. 

Droplet Tra jectory Calculation 

To calculate the droplet trajectory requires the 
numerical solution of eg. (4) . Eq . (4) is S^Olved in the 
cartesian coordinate system shown in figure 11. The T-y 
axis is used for' the trajectory calculation while the x^- 
y* system is usexl in the flowfield code. All inputs and 
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outputs to the trajectory code are in the x-y system. The 
initial conditions n eed e<! to solve eg. (4) are the droplet 
velocity and position in the free stre am. The particle is 
assumed to be travelling with the free stream at some 
finite distance in front of the airfoil, usually five 
chord lenqths. The initial y coordinate is selected so 
tne particle either strikes or misses the airfoil as 
desired . 

5q . (4) is a. second_or.der-, nonlinear, ordinary 
differential equation. Equations of this type are 
generally written in component form and reduced to first 
order for numerical solution. This results in a system of 
four simultaneous differential equations which can be 
solved by a step intoqration method. Kowever this system 
is stiff, and requires special numerical treatment for a 
stable solution. 

A stiff system has in its general solution eigenvalues 
which may be orders of magnitude different in absolute 
value and therefore each dominates the solution in 
different regions. If not handled projierly this leads to 
unstable solutions C&2]. This numerical formulation uses 
a variable step sixe, predictor-corrector scheme suitable 
for stiff systems by Gear [63,64], when compared with the 
Adams method on this system of equations, the stiff method 
reduces the computation time by at least a factor of two. 
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The system of differential equations can now be solved 
if a local velocity vector, u, and a droplet draq law are 
provided. The flowfield velocity calculation will be 
discussed in detail in tlhe next section. This program 
calls a subroutine which provides the velocity at any 
point (x,y) in the flowfield. Several droplet draq 
equations are available as discussed in Section III. This 
proqram uses the drag law of Langmuir f9] given in eg. 

( 8 ) . 

A trajectory calculation is termin ated when the 
particle* strikes the airfoil surface or misses and moves 
past the body . Polynomial fits to the trajectory and 
airfoil surface are used to determine the exact impact 
point, 6 as shown in fiaure 9, and the surface length S as 
iji fiqure 7. The tangent trajectories, figure 7, are 
calculated using an extrapolation procedure based on the 
impingement a-ngle. 

Usina this method the program can supply the Ay^ ,and y^ = 
Yq(S) needed to-caJ.. culate the local and overall collection 
efficiencies. All these calculations are controlled 
internally by the computer prog-ram, by error limits input 
by the user . 

Piqure 12 shows a typical yQ versus S plot generated 
by the program. The symbols are the results of actual 
dro plet traject ory calculations. These points are curve 
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fit usinq a cubic s-pline which forces the slope to zero at 
each ooint. This scheme for _spline fitting the Jq vs 

S curve must be modified for certain special cases. Pot 
large values of K the airfoil upper or lower surface, 
dependina on the angle of attack, may collect ice all the 
way to the trailing edge. In this case 8 does not equal 
zero at this limit of impingement and therefor the second 
derivative, rather, than the first, is set equal to zero at 
this endpoint. 

Another special case-results when an area of the 
airfoil, between the maximum limits of impingement, 
collects no ice. This results in a discontinuous Yq vs S 
curve. In this case the curve is fit in two pieces which 
are connected by a region of zero impingement, 8= 0. 

This second case occurs on airfoils with cusps, such as 
the MACA six series airfoils. Here the most forward 
region of the cusp may collect no ice tor large K’s and 
high a’s, while the aft segment does collect ice. This 
may also occur near the leading edge when time stepping 
leads to a concave region in the ice shape. 

The spline fit is then used to calculate the local 
impingement efficiency, 8» which is the slope of the 
curve, figure 13. The 8 distribution and airfoil geometry 
are stored on disc to be used for the ice shape 


calculation 
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Plowf ield 

Th^ velocities required for the solution of 

eq. (4) are generated using the Theodorsen method. A 
modified version of the flowfield code by Woan C60] is run 
once and the transformation results are stored on disc. . 
Input to the flowfield code are the airfoil coordinat es in 
the x’-y' coordinate system, figure 11. The droplet 
trajectory code reads in the results of the 
transformation. 

When the velocity at any x-y point is required, the 
velocity subroutine in the droplet trajectory code first 
must rotate to the x*-y* system, then transform the x*-y* 
point to the circle plane of the transformation. The 
transformation to the circle plane is nonlinear and 
therefore a Newton-Raphson iterative technique is- used. 
Once in the circle plane the velocity calculation is 
straiahtforward . Note that this method calculates the 
velocity from the transformation at each point required by 
the step integration differential equation solver; a 
matrix of stored velocities with an interpolation scheme 
is not used by this program, 

A second program by Moan [f>0} is available to 
calculate an inviscid C Cm, and if desired. This 

program also generates a Cp plot which is useful in 
ensuring smooth airfoil coordinates. 
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Ice Shape Calculation 

Eq . (36) must be solved for ^ to deterndne the rime 
ice shape. The ice shape prediction code reads in B as a 
function of S, 6 (see fiqure 9) as a function of S, and 
the the airfoil coordinates f rjom_Lhe_disc file written by 
the droplet trajectory code. The accumulation parameter, 
Ac, is the only physical variable read in directly by the 
proqram. Internally the pr^ram must calculate <t> , ’ll , 
and e, fiqure 9, and either the surface radius of 
curvature, r, or the effective radius of curvature, r», in 
order to solve for I an d -calculate the ice shape 
coordinates , 

For normal ice growth, eg, (37), the surface radius of 
curvature, and the direction of the outer normal, e, are 
needed as a function of S, Both terms can be found from a 
polynomial fit to the airfoil coordinates. For airfoils 
with roudh coordinates, e is calculated at the droplet 
impact points and e vs. S is fit using a cubic spline. 

From the cubic spline e and r ( r = ~d£/de ) can -be 
calculated 'at any S location. This procedure provides 
smoother values ‘of e and r. 

For non-normal ice growth, eq . (36), r», and il; must 

be determined. For the tangent case 6 is known at each 
particle impact point and £ can be calculated f.ro» a 
polynomial fit of the airfoil. Then (ji ( = e ♦ n/2 - 6 ) 
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versus s can he spline fit and r», eg. (38), and 41 can 
be found for any S location. The code allows for ice 
growth directions-Othe r th an normal, e, and tangent, ip. 

By redefining the angles tp and 4> to be measured with 
respect to the assumed ice growth direction, instead of the 
trajectory tangent, the same method that was used for the 
tangent cajie can be used here. Then the assumed ice 
growth direction 'P can be chosen arbitrarily. 

With S, and the direction of growth determined, each 
airfoil coordinate affected by -the ice is recalculated. 
This generates the iced airfoil coordinates. A 
trapezoidal integration is used to determine the ice shape 
area to be checked aaainst the exact area, eg. (35) . The 
original and iced airfoil coordinates are written on disc 
for input to the next code. 

Iteration and Smoothing 

Airfoil analysis codes are in general very sensitive 
to-b^^ the first and second derivatives of the airfoil 
shape as provided by the input coordinates. After the 
airfoil has been iced, the coordinates are often too rough 
to run well in these programs. Existing airfoil 
coordinate smoothing programs were not designed for, and 
therefore can not handle the type.s of airfoil shapes that 
result from the icing analysis. Therefore, a coordinate 
smoothing program was written specifically for the icing 
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problem . 

The smoothing is accomplished by force fitting a 

polynomial of the form [65} 

^ + • • • • + C^X + Cq (39) 

to both the upper and lower surface of the ice shape. The 
exponent p is a fraction to allow the notching of the 
leading edge radius of the ice, and N is the order of the 
polynomial. The desired first derivative is automatically 
satisfied at the leading edge, and the function is forced 
to match the slope of the airfoil surface just, aft of the 

f) ice accretion. 

h 

i An additional smoothing routine is available when the 

i 

( ice shape is not of the form of eq . (39} . In this case 

the ice shape is essentially smoothed by hand with the 
help of an interactive computer graphics program. The 
program displays the original airfoil leading edge and the 
new iced airfoil shape. By using the_ cursors the iced 
airfoil coordinates can be adjusted to provide the desired' 
smoothing and coordinate distribution. 

When ti me -stepping an ice build-up the smoothing 
program is available to generate iced airfoil coordinates 
to be used in the flowfield code. Depending on the value 
of the accumulation parameter, coordinate smoothing may or 
may -not be required for every time step.. On the last time 
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V. AEFODYNftMIC AJUILYS15 

Th^ roost, serious of icf^ foirro3t.ions on airfoils 

aro Iho roaurlions in roaxirouro lift coefficient and a 
sicrnificant in draq. JlijDe changes the airfoil 

q<x)metry and adds roughness to the airfoil# These two 
effects are primarily responsible for the change in 
airfoil perforroanco due to riroo ice. Existing airfoil 
analysis codes are able to anuly7-e the iced airfoil shape/ 
but do not properly handle the roughness effects# As a 
result/ the effect of the change in airfoil shape and 
surface roughness roust be handled separately# The new 
airfoil shdw will be handled analytically/ while the 
rouqhnesf^ effects will he accounted for using empirically 
based corrections • 

Te r Shape Analysis 

Riroe ice accretions are streamlined in shape but do 
not blend smoo*thly into the airfoil shape. In addition 
the shape itself may not be ^smooth** with respect to the 
requirements for good airfoil leading edge geometries. 

Due to the qooinetry of the ice shape severe adverse 
pi'TSsure gradients occur in the leading, edge region. 

These gradients trigger the. -f.oxma.tion of jarnall zones of 
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separa-teil flow (separation bubbles) which at higher angles 
ol attack may lead to massive separation and stall. While 
-saarface roughness may ul; o triggger premature stall. This 
analysis assumes that the reduction in maximum lift 
coefficient of iced airfoils is due to the change in 
leading edge shape alone. 

The Kppler [61] airfoil analysis code is used to 
predict the effect of the ice shape. The code uses a 
sophisticated potential flowfield model of distributed 
surface singularities with parabolic strengths on~curved 
surface panels. The version of the code used has been 
modified to include the compressibility effects' on the 
potential flow. Under this potential flow an integral 
boundary layer method is used to calculate the skin 
friction. A louuhness factor is included but its only 
effect is to cause early transition from laminar to 
turbulent flow. K special feature of the program is an 
approximate calculation of the maximum lift coefficient. 
The lift is ca-lculatod by usina the two dimensional lift 
curve slope and a corrected absolute angle -of attack. The 
correction reduces the angle of attack based on the size 
of the separated zone- 

The airfoil analysis program is used to analyze both 
the original airfoil and the airfoil with the rime ice 
shape-. The program provides the lift, drag, and pitching 
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moJiM^nt copf ficif'Hts for both casps as well as detailed 
pressure distributions. The draq prediction for the ice 
P-haoe must still — lie corrected for roughness effects. 


P.ouahnesK- Effects 


— roughm'ss caused by rime icing is large compared 
to the boundary layo'r thickjjess. This roughnoss not only 
increases the local skin friction, but it can remove a 
censiderahle amount of kinetic energy from the boundary 
layer. This increases the skin friction drag and adds 
pressure drag -due to the base drag of the rouahness 
elements and the reduction in pressure recovery due to the 
thic)teninq of the boundary layer [2]. This reduction in 
pressure recovery can lead to premature stall due to 
boundary layer separation at lower than expected angles of 
attack 

In a recent paper Drumby [66] has compiled the 
existing data on the effect of roughness on maximum lift 
coefficient. This summary is shown in figure 14. The 
data shows the rather dramatic reduction in maximum lift 
due to relatively moderate levels of roughness. Also 
presented in Brumby's paper is u good discussion of the 
ojerational aspects of wing surface roughness. Although 
figure 13 will not be used directly in this analysis, it 
does provide a good check On the analytical results. 


Grav [31 pr(^s«^nted an empiricf.l corrf*lation to pr<*dict 


dratj incr^’Hien tc dm-? to airfoil icinp. 

4Cd -[8.7 X 10-5^,J[;;;;;(32-t)“'^](i + 6{(1 + 

2.52r*^‘^ sin^ 12 a-) sin^|543V^^ — (40) 

This equation was, however, developed primarily for the 
qln 7 .e ice case which was felt to be the more serious 
problem. Tne rorrelation is linear with time which .does 
not accuratelv represent the time data. Therefore a new 
correlation is neenoo which is developed specifically for 
the rime ice case. 

Ttip amount of oocd data available for the draa of 
airfoils with rime ice is very limited. Therefore the 
problv^m was formulated to ta)te advantaqe of the data on 
airfoils with leadinq edqe roughness. (When good iced 
airfoil drau data is available, this correlation could be 
easily modified to inclnue this new information.) Figure 
If) shows the (!ra<i increase vexjsaxs ire accumulation (a 
function of time) for both daze and rime conditions [36]. 
Note that the increase in draa for the qla.ie ca.se. i s 
aporoximatelv linear as Predicted by Gray's eq. (40). 
However for the rime case, the dran increases rapidly at 
first, then levels off and increases linearly a t -a reduced 
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ratp. This analysis, as shown by the dotted line, ignores 
the Jjiitial rapid increase and matches the linear section 
assuming a step increase in drag as soon as icing begins. 

The Latenoapt of the linear drag law proposed can be 
obtained from figure 16- These empirical curves were 
obtained from published experimental results on airfoils 
with leading edge roughness. Note that different types or 

families of airfoils are affected differently by _ 

roughness. These differences are due primarily to the 
amount of laminar flow the clean airfoil experiences. 

Gray allowed for this change by including terms based on 
the airfoil leadino edge radius. Given a particular 
airfoil, figure 16 can be used to estimate the step drag 
riJ5o. A value of K/c * 0.001 is representative of the 
initial roughness of the ice. 

With the constant term in the proposed di’ag 
correlation determined, the form of the time dependent 
term must be developed. The independent variable must be 
dimensionless to remove the scale effect. For example, 
two airfoils of different chord lenaths which have the 
same scaled ice accumulations should have the same 
increase in drag. Represontina the drag rise as a 
function of ice accumulation would however give these two 
air foils- different drag increments. A better choice of 
the independeiut- variable is the dimensionless collection 
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parameter, AcE. This is just the cross sectional area of 
the ice shape divided by the projected height of the 
airfoil. Here the initial value of E from the theoretical 
analysis is used and note Ac is linear with time. 

Figure 17 shows some of th e callable rime ice data 
plotted versus the collection parameter, AcE. Note that 
for all the airfoils the slope of the curve is the same. 
The predicted results shown on figure 17 u-se the values . . 
from figure 16 for t-he AcC = 0 drag increments. 

Exoressing the res-ults. of figures 16 and 17 in_ equation 
form 

LC^ = .0l[l5.80 ln(^)+ 28000 A(,E + l] (*♦!) 

where I is the constant which depends on the airfoil type. 
Table 3. 


Table 3 Constants For The Drag Equation 


Airfoil Type 

Drag Constant, I 

Typical k/c 

4 and 5 Dioxt 

184 

.001 

66 Series 

218 

.001 

Series 

232 

.001 - 

65 Series 

252 

.001 

66 Series 

290 

.001 


The new drag of the iced airfoil is then given by 


Note that in all cases ACd is based on the G»d for the 
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hydraulically smooth airfoil at the given angle of attack. 
This removes a possible source of error since all models 
Biay have different roughness levels due to the 
construction techniques or condition of the surface. 
Analysis Procedure 

Tiie aerodynamic analysis can be suminarizeU as: 

1) Calculate the icing characteristics and rime ice 

shape usinq the procedures described in Sections 
111 and IV . . 

2) Dse the airfeil code to analyze the clean 
airfoil 

3) use the airfoil c-ode to analyze the smooth iced 
airfoil to predict the change in maxiaium lift 
coef ficient 

4) Use eg. (41) to correct the drag analysis for 
roughness effects 

Step 1 not only predicts the ice shape but the 
collection efficiency, E, which is needed to determine the 
drag in step 4. Next the clean airfoil performance is 
analyzed to provide a baseline and also to generate the 
value of Cd which the correlation of step A is based. The 
smooth ice shape is then analyzed using the- airfoil code 
to determine the maximum lift and pitching moment. 

Finally the empirical corrections are made to yield the 
effect oh drag due to the rime ice. This correlation 
relies upon published data and the results of steps 1 and 
2. This method for analyzing the aerodynamic effects of 
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rime ice on airfoils uses the analytical methods which are 
available or have been developed here, and supplements 
these with empirica I* resu ljts_Sfhen needed. 
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VI RESOLTS AND DISCUSSION 


Thf' purpose of this study was—to develop an analytical 
method to ajialyze the rjime icinq~of airfoils. Therefore, 
this section deals primarily with the validation of. this 
methoo . The analysis will be. coinpar.ed to other analytical 
results and to the experimental data which are available 
or were generated specifically for this validation, la 
addition, limited use of. the method has been made to 
analyze the effects of certain parameters on icing rates. 
Trajectory Analysis Validation 

Langmuir [9] first formulated the droplet trajectory 
equation for numerical solution on a differential 
analyser. Several calculations were made for the case of 
a circular cylinder, sinc€* this flowfield can be expressed 
in closed form. Lanamuir's results wt-re often used as 
test cases for the NACA and other trajectory calculation 
methods . 

Table U is a summary of some o4 the analytical 
predictions of icina rates on circular cylinders. 


Table 4 


m 

m 



e>3- 


Comparison of The Present Method to That of 
Langmuir [9] and Lozowski [67] for Cylinder 
Icing Pates 



Langmuir 

Lozowski 

Bragg 

Case No. 

1 

2 

1 

2 

1 

2 

R 

600 

100 

60 0 

100 

600 

100 

K 

18 

0.5 

18 

0.5 

18 

0.5 

Vx 

1.056 

0.494 

1.056 

0.477 

1.026 

0.425 

Vy 

0.193 

0.725 

0. 196 

0.650 

0.19b 

0.623 

E 

0.819 

0.156 

0.814 

0.170 

0.812 

0.155 


0.885 

0.34 8 

0.898 

0.376 

0.900 

0.363 

•^tnax 

®m 

79.8 

34.2 

79.5 

35.6 

79.1 

34.it 


Included in the table are the results of Langmuir and 
Blodgett [9], Lozowski and Oleskiw [67], and the present 
method. The results of two test cases are shown. Here Vx 
and Vy are the dimensionless velocity components of the 
tangent trajectory particle as it strikes the cylinder. Bmax 
is the maximum impingemont efficiency which for a circular 
cylinder with no circulation occurs at 0=0 degrees. 

Here 0 is the angle which defines a point, on the cylinder 
with 0=0 being the forward stagnation point. The limit 
of impingement for the symmetric case is then 0jjj7 

As indicated, all three methods agree very well on the 
first test case, table 4. The agr€;ement is within, one 
percent on th(' value of E and 0j^ , while the values of Bmax 
are within two percent. However for case two, while the 
agreement is good, there are some more significant 
differences. Case 2 is a more severe test than Case 1 
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since the value of K is_almost 20 tines smaller. This 
small value of K results in particles which are much more 
affected by the flowfieid and therefore their trajectories 
are more difficult to calculate accurately. Fere 
Langmuir’s method and the present method agree closely, 
while LozowsHi’s calculations are about ten percent higher 
in collection efficiency, 

^*he source of the differences is not obvious. All 
three methods use different equation solvers, drag laws, 
and flowfield models. The allowable errors in the 
numerical schemes may also be different. The present 
method was run with the allowable error in E not to exceed 
one percent. No error tolerances were reported by 
Lanqmuir or Losowski. The most likely explanation of the 
difference in Lozowski's calculations is the drag law 
chosen. Since these particles do have low inertia, a 
small change in the assumed Cd could have a large effect 
on the results. 

Droplet trajectory calculations can also be compared 
to early NACA results for impingement on a NACA 65A004 
airfoil. Figure 18 shows the early calculations [19} 
compared to the present method for the airfoil at zero 
degrees angle of attack. The comparison is quite good 
considering the errors involved in the early caJ.culations . 
Brun [ 17 ] estimates the error in 6 for "the NACA method to 


about ten pprcent. This is clue to the se?ere velocity 
qradients around the small leadiny cdqe radius and the 
difficulty in curve fitting, and dettirJLining the slope of, 
the Yq versus S curve to get 3« The present method 
performs this calculation routinely to within one or two 
percent . 

Pecently analytical results of airfoil droplet 
impingement have been published by Lozowski and Oleskiw 
[ h7 3 , Lozowski *s general numerical, scheme is the same. as 
the present analysis, while the details of the solution 
varies in several areas. Figure 19 is a comparison of 
Lozowski *s results and the present method for a NACa 0015 
airfoil at eight doareos angle of attack. The results of 
the two method? arc In good agrcjcmont in all areas. The 
limits of impingement, , and the 6 curve itself are 

practically identical. Lozowski's reported collection 
efficiency of 0.501 seems high when compared to the two 
curves and the value of 0.473 for the present method. 

Figure 20 shows a similar comparison at a slightly 
different condition. However here Lozowski [67] has 
included the Eassett unsteady memory term which was 
dropped from the differential equation used in this 
method. The comparison is still good, with Loz< wbki's 
results showim --re droplet impingement. The addition 
greatly complicates the droplet traiectory calculation and 
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rpRiilts in only a small chanqo in g. This correction is, 
however, within the the error caused by the difficulty in 
measurinq the droplet size distribution in a cloud, and 
also the error inherent in a sphere draq curve fit. 

Limited experiaiental data is available for water 
droplet impinyement rates on airfoils [30]. These data 
were taken usinq the d-ye tracer technique in the NACA 
Icdnq Research Tunnel. Impingement data taken on a HACA 
6h-212 airfoil at four degrees angle of attack are 
compared to the theoretical results of this method in 
fiqure 21. The comparison between the theoretical and 
experimental results is quite good. The absolute value of 6 
from the experiment may not be accurate due to the 
problems in the calibration of the free stream conditions. 
However the limits of impingement and overall character of 
the curves compare very well. It should be noted that in 
the experiment the droplets were not of a single uniform 
size as was assumed in the present calculation. This 
point will be discussed in the next serrLion . 

The present method and computer code for calculating 
droplet trajectories and ultimately impingement rates has 
been compared to earlier works. Results from two very 
early analytical methods and a recent Canadian method 
compare very well to the present results. These 
comparisons were made on both airfoils and cylinders. 
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Comparison of tht- present, m ethod to experimental results 
was also shown to be very qood. The pjresent method has 
therefore been shown to be valid and yield ver y ac curate 
droplet impinqement results, 

VMD Approximation 

Actual icina clouds contain a distribution of water 
droplet sires. Fiqure 22 shows th<^ resultinq 6 curves for 
droplets from 10 to 50 microns impina inQ on a NACA 0012 
airfoil at an anqle of attack ©f five deqreep. The 
traiectories of the smaller particles are dominated by the 
draq term in the di'^ferential equation since the inertia 
is small. The droplets fjollow the streamlines more, 
closely and therefore few imninae on the leading edge. 

For the la-r.q-Gr droplets tne inertia term dominates and a 
large percentaae of the particles impinae on the airfoil 
leading edge. Note that the area under the B curve is 
proportional to the total mass striking the airfoil, 
therefore clouds of larger particles will increase the 
mass of ice accreted . 

nsino the method of Section TIT and a Langmuir r 
distribution of particle sizes, a $ curve fo.; the entire 
cloud of nonuniform droplet sizes can be predicted, figure 
23. Also depicted in figure 23 is the B curve for a single 
droplet size, the volume median diameter, VHD. The ?PID'is 
the droplet diameter for which half the volume of water in 
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thf> cloud is made un of drop^lots larger than the VHD, and 
half the volume from droplets smaller than the VHP * As 
seen in the figure, the VMD 3 curve is a very good 
approximation to the actual icing cloud results. The VMD 
anproxiination slightly over predicts the 3 max and has 
reduced majcxniuni limits of impingement. However these 
errors are acceptable in exchange for the reduction in 
computer time. lonoring the droplet size distribution 
effects saves an order of magnitude in computer time by 
reducing the number of droplet diameters which must be 
run. In avddition it eliminates completely the 
calculations needed to combine this information into one 
r-ijrvp. Therefore, unless stated otherwise, all 
impinqement calculations presented here will use the VMD 
anofoximation . 

Scalina Parameter Validation 

The simplified similarity parameters and R have 
ht'en derived in Section III. Both parameters combine Ry 
and K into a sinole dimensionless quantity which greatly 
simplifies the icina problem. Those parameters can be 
used to facilitate data presentation and to define test 
ronditionr. for scale model tests. Here expecimental and 
numerical data are used to compare and evaluate the 
modified inertia parameter, K^, and the trajectory 
nimilnritv parameter, K:._ 


w J' |t»' • -I t. •■' «' ' *' ■■• 11 ^ 
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Historically icing data has been presented using the 
modified inertia parameter. The deqreo to which 
Compresses this data to a single curve provides a measure 
of the accuraci_of the approximation. Piqure 29 shows the 
airfoil collection effici»>ncy, t, for three different free 
st-reani Re-ynolds numbers and for various values of K 
plotted versus Kq. The results are from an early NACA 
analytical study [17] of a NACA 65A004 airfoil at fjour 
deg reps_ajuale of attack. The same data are plotted as a 
function of K in fiaure 25. Here C is taken as one and 
Y = 0.35 as discussed earlier-. Both parameters reduce the 
data toward a sinqle curve, but the K parameter shows 
somfMfhat less deviation from the curve. It is not clear 
from these results if the scatter in the data is caused by 
the- similarity oarameter approxima-tion, or if the error is 
in the numerical results for E. 

To attempt to resolve this uncertainty the present 
(Jroplet trajectory code was used to generate similar data. 
Htto a NACA 0012 airfoil at zero degrees angle of attack 
was analyzed at three different values of Py and five 
values of K. These results, plotted as a functioa of Kq 
and K, are <jiven in figures 26 and 27, respectively. Here 
both and K do an excellent job of reducing the data to 
a sinqle curve. This suggests that the scatter in figures 
29 and 25 is error in the early numerical data, and not a 
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reflection upon the accuracy of ani3 R. 

Both the modifiefi inertia parameter and the trajectory 
similarity oarameter simplify the droplet trajectory data 
presentation. An additional numerical check on the 
validity of the parameters can be made by comparing scaled 
d-rople.t impinuGTnent efficiency curves-. The results of 
using Kq and R as scaling parameters for a one-sixth scale 
model are shown in figure 20. These curves were generated 
using the method and computer code described earlier. 

For scalina droplet trajectories the K' parameter hajs a 
definite edge over since y may be optimi2ed for each 
droplet size (the VMD If a distribution is considered). 

The procedure used for determing Y described earlier- 
yields a Y of 0.30 for the 15 micron full scale droplet 
and 0.39 for the 30 micron size droplet size droplet. The 
values of R and K used as well as the droplet diameter, 
are given in table 5. 
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Table 5 Scaled Variablep for Analytical Icing Test 
Usiiug— K. and K 



Full Scale 

One-Sixth 

Scale node! 



K 

Ko 

6 ( ym) . 
K - 

15.0 

115.6 

0.0393 

5.23 

40.30 

0.0286 

5.05 

38.93 

0.0267 

6(ym) 

30 .0 
231.2 
0.1572 

9.86 

75.97 

0.1018 

9.6 0 
73.98 
0.0966 


Note that fon this example it «as assumed thaJ:;._only 
the particle diameter would be changed to provide the 
scaling. All other variables such as the aircraft 
velocity, droplet density, air density, etc would be held 
constant. This yields an equation for the droplet 
diameter of -i- 



The important results of the scaling comparison of figure 
2B are summar izeja_JLn. tjauble. 6 -. 


abl<= 6 R(-sults of Ihe Droplet Trajectory Scaling 
Comparison 


full Scale 


One-Sixth Scale 


V- 

^max 


6* 1b pDl 


0.05L)b 

0.332 


0.0557 

5.331 


6* 30 yin 


0.050S 
0.32 3 


0.173 

0.560 


0.174 

0.563 


0.166 
0.56 3 


While K docs a reasonable job of reproducing the full 
0 

scale ts^ jectories, the added flexibility in the K 
paraiiiGter allows for an excellent trajectory scalinq. No 
experimental results are available to evaluate the 
similarity iaramotcrs for the airfoil icing scalino 
problem .. 

However, recently published experimental results by 
Orinsbee and Fragg [54] are available for a similar droplet 
trajectory case. In these tests conducted in the NASA 
Lanoley Vortex Research Facility, three geometrically 
scaled aoricultural aircraft models were used to inject 
scaled spherical particles into the model wake. Dsinq th» 
com.plGte set of similarity parameters for the droplet 
dynamics Rpj, K, and Fr results in a unique scaled test 
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particle o£ low density and large diameter. Relaxing the 
constraints on the scaled particles by repl acing Ey and K 
by K yields an infinite number of candidate test — 

particles. This greatly simplifies the task of obtaininq 
tfie test particles. - While Ormsbee and Bragg did not use 
in the same form as it was derived here, their method 'is 
completely equivalent in that they made a similar scaling 
aporoximation . 

In these tests a hy {othet icai full scale aircraft and 
droplet test conditions were- chosen. These were then 
scaled to determine the equivalent test conditions for a 
O.fO, 0.15, and 0.20 scale model. Table 7 shows the full 
scale and model test conditions while the particle 
trajectory results are summarized in figure 29. 


Table 7 Scaled Physical Variables for Droplets 
in Aircraft Wake 




Model 

Scale 



0.10 

0.15 

0.20 

1.0 

Rinq semi span, :n 

1-.22 

1.83 

2.44 

40.0 

Model velocity, m/sec 

16. ft 

20.6 

23. ft 

53.3 

Altitude, m 

.622 

.933 

1.24 

20.4 

Angle of attack, deg 

2.00 

2.00 

2.00 

2.00 

Particle diameter, ym 

105 . 

125. 

105. 

490. 

Particle density, g/cm^ 

2 .42 

2.4 2 

3.99 

1.00 


Presented in the fioure is the lateral transport of the 
particles bv the wake vortex system as a function of the 


lit: 
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initial injector location. For all three models the 
lateral tra nsport-oi-^he scaled particles is the same, 
verifying the K scaling analysis. Scaling_JLejsdLs.we.re also 
conducted [ 54 ] in which other lift coefficien ts , aircraft 
altitudes, and full scale droplet sizes were used and in 
all cases the particle trajectories scaled well. 

Trajectory Results 

Although the objective of this study was to generate 
rime ice shapes and evaluate their aerodynamic 
performance, the trajectory calculation' ilone provide 
much useful information. The droplet trajectory computer' 
program can be used to conduct a sensitivity analysis and 
provide physical insight into the impingement process. 

The information provided by the analysis such as the 
overall collection efficiency and maximum limits of 
imninoement can be used directly in the design of ice 
protection system.s. 

Fioure 30 shows the paths of water droplets around a 
NACA 0012 airfoil. Trajectories are shown at both zero 
and five degrees angle of attack. Note that at five 
degrees the droplets which impingement on the airfoil 
start out below the airfoil in the free stream. This is 
of course due to the upwash in front of a lifting airfoil. 
The Darticles which miss the airfoil by pass-ing over the 
leading edge gain a large amount of Kinetic energy in the 
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leading edge region. These particles are therefore less 
influenced by the flowfield over the aft part of the 
airfoil. Although little quantitative information is 
obtained froa— the trajectory plots, some physical feel for 
the problem can be gained from them. For example, droplet 
trajectory plots proved very valuable in identifying _ 
regions on— t.hjt?-_ai.r.f.Q.H_ where no particles hit the surface. 
This led to modifications in the spline fitt ing program as 
described' in Section IV. 

The effect of airfoil anqle of attack on droplet 
impingement efficiency is shown in figure 31. As expected 
the area of impingement moves more toward the lower 
surface as tl '? angle of attack is increased. Also the 
area under the B curve, the total mass collected, 
increases with angle of attack. A slight change in the 
location and value of B max, the maximum local 
impingement, occurs with the increase in angle of attack. 
This effects the sh-ayie of the leading edge ice shape which 
may cause large differences in the aerodynamic performance 
of airfoils iced at different angles of attack. 

A sensitivity analysis may also be pertomied by 
varying the value of K, the inertia parameter. Varying K 
while holding Ry constant corresponds physically to 
subjecting airfoils of different .chord lengths to the same 
icing conditions. Note that the airfoil chord, c, appears 
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in the denominator of K , so reducing c increases the value 
of K. Increasing K while holding Ry constant means that K 
increases linearly with K. 

Figure 32 demonstrates the effect of varying K, or 
equivalently R. Here the case of a NACA 0012 of chord 
six, three, and two feet (increasing K) is shown. Then as 
the airfoil chord decreases the overall -collect ion 
efficiency, area under the 6 curve, increases. Since the 
smaller airfoils have more severe velocity gradients near 
the leading edge, the droplets are not able to follow Hue 
streamlines as well, and more droplets impinge on the 
airfoil. This is observed in flight when tail surfaces, 
because of their smaller chord, accrete proportionately 
more ice than the main wing. It is interesting to note 
tiiat for the range of K represented by figure 32, K = .008 
to .025, the collection efficiency as given in figure 27 
is almost linear. In fact for this special case as K 
increased 200 percent, su did the collection efficiency, 

E. 

Figure 27 also represents anothsr use of the method. 
Osina K as the independent variable, the initial icing 
rates anJ other results may be generated to evaluate the 
susceptibility to icing of a particular airfoil.. Here 
only E is nresented as a- function of K, but a complete 
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airfoil analysis would include plots of S 

and the actual 3 curves, _ “ - JT ♦ 

Ice Shape Calculation 

Before the aerod ynainic tjerfurmance of an airfoil with 
rime ice can he determined^ the ice shape must be 
accurately predicted . TJiis involves the time-ste pulnu 
procedure outliiiod in Section III, Havino shown JtlraJLJt he . 
initial icinq rates predicted by th^ method are valid# the 
accuracy o-f the time^stej ping model to predict rime shapes 
will now be examined. 

First the assumed direction of growth out from the 
airfoil surface must be determined. Figure 33 shows a 
normal and tangent growth predicted from the same initial 
droplet impingement information. Both shapes represent 
one large icing step, that is no time-stepping was 
performed. The predicted tanaent shap^» grows out into the 
oncoming droplets. With its increased maximum growth and 
reduced leadina edge radius it 4ias the general shape of a 
measured ice accretion. However physical intiiition would 
suggest the normal growth to be th^ correct mode- In the 
limit as the icing time goes to zero, tne tangent growth 
approaches the same shape as the normal mode. It is felt 
that the normal growth model is the physically correct 
solution for a time-stepping procedure. The tangent 
growth appears to be an approximation to the time-stopping 
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niHthoil — vil] bp more obvious later* 

timo-steppinq proce(^ure is den onst rated in figures 

34 and 35 on a modified NACA 64-215 airfoil at a cruise 
condition. Here the angle of attack is 0,3 degrees and 
= 115.6 while K = 0.044, Three time steps-were taken, 
each representing f4.ve minutes- of icing with the 
accumulation parameter, Ac, equal to 0.0133. Figure 34 
shows the predicted ice -shapes— from the 6 curves of figure 
35. Note that in the time-stepping method first the 
impingement efficiency is calculated on the clean airfoil, 
step 1 figure 35. Then this g curve is used to predict 
the first ice shape figure 34. Thp flowfield is then 
recalculated, the step 2 6 curve generated, the new ice 
shape 2 predicted, alid the iteration is continued. 
Therefore figures 34 and 35 are intimately related. 

Examing figure 35 the changing airfoil shape is seen 
to have a significant effect on the droplet impingement 
characteristics. This supports the need for a time- 
stepping approach. The change in the impingement values 
with time are summarized in Table 8. 
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Tablo 8 Time Step Parameters 


Step 

Ac 


^^o 

AS 

1 

0.0133 

0.358 

0.00983 

0.04 95 

2 

0.0133 

0.411 

0.00909 

0.0379 

3 

0.0133 

0.472 

0.00910 

0.0342 


The maximum inpinperacnt efficiency 3 max, increases with 
time while the iced surface length on the airfoil, AS, 
decreases. The_oJBLe.rall collection efficiency decreases 
slightly for this case. Another interestin g feature is 
the development of the second peak in the curve on the 
third time step. 

All these effects of time are also reflected in the 
predicted ice growth, figure 34. The increase in 3 max 
and redaction in As generates the reduced leading edge 
radius of the ict' and the imore pointed shape. The second 
peak in the 3 curve results in the reflexed upper surface 
"bump” on the third time step- The effect of time- 
stepDino is then apparent from the change in the curves 
frcm steps 1 to 3. 

The accuracy of the time stepping model will be-a 
function of the size of the time step taken. Figure 36 
shows the predicted ice shape for the same arcndifled 64-215 
airfoil with one, three, and six time steps. Here the 
correspondinq Ac's are 0.04, 0.0133, and 0.067 
respectively. The 3 curves for the six time step-case are 
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qiven in figure 37. A significant change in shape is seen 
between one am^hree steps, while the change front three 
to six steps is relatively smaJl. In fact the change in 
shape froJD. three to six steps is probably due as much front 
nutnerisa.1 error as front an improveniGnt in the physical 
modeling . 

A s-imilar study on the ef£ect.jof step size was 
conducted using a NACA 65A413 airfoil. The airfoil was 
analyzed at one degree angle of attack with Ry = 147 and 
K = 0.118. The length of the icing, ejicounter is eight 
minutes, which for the free sireaiii conditions assumed, 
gives an accumulation parameter for the total time of 
0.044. The predicted ice shapes for one, two, and four 
time steps are shown in figure 38« Figures 3? and 40 are 
the corresponding 8 curves. Here the shapes do change 
from the two to four time step case. The maximum amount 
of ice growth, and the shape of the ice near the limits of 
impingement, do not agree. The four time step case has 
essentialXy. taken mass from near the limits of impingement 
and shifted'it forward by extending the leading edge 
growth . 

From these two cases, and other experience with the 
method, some guidelies in selecting the step size can be 
formulated. The critical area is the region near the 
maximum limits of impingement. The ice’ in this region 
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should not be allowed to grow more than that which 
generates a shape that blends in smoothly with the 
airfoil. A rule of thumb is that the maximum growth in a 
single step should not exceed one-half of one percent 
chord, X = 0.005. This corresponds roughly to holding 
< 0.005. The allowable step size is actually a 
function of the leading edge geometry and the shape of the 8 
curve. With airfoils with small leading edge radii 
requiring the smaller stop size. The rule of thumb givjan, 
however, provides guidance in selecting nn acceptable step 
size. The lower bound on the step size is governed by the 
amount of computer time required per step and the 
accumulation of numerical error. Error accumulates 
primarily due to the coordinate smoothing process. The 
smoothing required is due in part to the discontinuous 
surface radius of curvature ol some airfoils and this 
problem is aggravated as more steps are taken. From 
experience the step size suggested appears to be an 
optimum for reducing computation time and increasing 
accuracy . 

With the time -step ping procedure established, this 
method for predictina ice shapes can now be compared to 
some experimental results. Experimental tests completed 
recently in the NBSA Icing Research Tunnel have generated 
experimental rime ice shapes for the modified NACA 6<4-215 
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airfoil [68], The ex peri mental rime ice shape for the 
cruise condition is compared to the present analytical 
method in figure 4JL. Tie ice accretion is small and 
therefore only the first one percent of the airfoil is 
shown, the experimental _shape and the time stepped 
preddction (from figure ifi.l_acuii,pare very well." The no 
time step case is also— shown to demonstrate the 
improvement in the prediction- when the time effects are 
included . 

The results of this comparison, and other test cases, 
permit some important conclusions to be drawji— The time- 
stepping was done fOr this case assuming normal ice 
orowth. Comp urinq figure 4 1 to the normal and tangent 
growth in figure 33 a similarity is seen. The tangent 
growth has the same general shape as the time stepped 
prediction. This suggests that the tangent growth is an 
approximation to the time effects. Also note that the 
time stepped shape pre4icts the reflexed upper-surface 
reoion and bumf)- as seen in the experimental shape. The 
upner node resembles the beginning of a second horn as on 
a glaze ice shape. However here it occurs solely as a 
result of the time effects on the f lowfield-and droplet 
dynamics. While glaze ice growth is certainly a 
thermodynamic orocess, this result suggests that 
impinuement characteristics may also be very, important. 
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For accurate alaze ice shape predictions, the time effects 
on the impingement rates should also be considered* 

This method has also been compared to the experimental 
results reported by Gray [37] on a NACA 65A00h airfoil. 

The airfoil is at two deorees angle of attack, R » 113 
and K = 0.J41, and the icing time is five minutes, Ac = 
0.0215. The experimental and analytical ice shape is 
shown in figure 42. The time -step ping improves the ioe 
shape prediction over the no time stepped case, but the 
shape is off considerably along the lower surface. 

The overall collection ef ficiency. "parameters, however, 
compare very well. Table 9. 


Table 9 Comparison of Theory and Experiment 
on the NACA 65A004 Airfoil 





Experiment 

Analysis 

W (Ih ice 

/ ft 

span) 

0.404 

0.331 

F.— - - 




0.208 

0.162 




0.0035 

0.0040 

St 



0.090 

0. 10 

Ice Area, 

ft^ 


0.0207 

0.0062 

Ioe Density {% 

o 

CM 

a: 

0 

31. 

85. 


The total masr collected, the collection efficiency, and 
the limits of impingement are very close for both the 
experiment anvi the analysis. However large d iscrepencies 
occur in the cross sectional area of the ice and the ice 
density. The error arises due to the assumed ice density. 
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8b percent the density ot water. This value is within the 
EAnge of 75 to 91 reported by Wilder [41] and close to the 
value of 89 used by LozowsKi and OlesKiw [67]. All these 
values are far from the 31 percent measured in the icing 
tunnel test. 

The very low measured value of ice density can be_ - 
attributed to the formation of rime feathers on the lower 
surface ice shape. These ice formations can be seen in 
the photographs and sketchs of reference 37. Rime 
feathers are thin layers of ice separated by layers ot air 
which sometimes f-orm during rime ice accretions, the 
occurence of rime feathers, which drastically reduces the 
overall ice density, is difficult to predict. These 
feathers cause the elfe»7tive ice density to be a function 
of S. If the correct ice density could have been used in 
the prediction of figure 42 the agreement would have been 
much betrter . Tbe prese nt method does not handle the rine 
feather case. However, when methods are available to 
predict the formation of rime feathers, this could easily 
be incorporated in the procedure. 

Aerodynamic Analysis 

The details of the prediction of iced airfoil 
performance is given in Section V. As noted there, little 
aerodynamic data is available for use in verifying the 
method. Therefore an airfoil test was performed on a 
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simulated rime ice shape to generate data for this 
purpose. the analytical method will first be compared to 
the simulated ice shape data. Then for the ioe shape 
- ■ predictions already distnjssed , the predicted airfoil 
performance will be compared to the experimental data. 

The simulated ice shape test was conducted on a NACA 
65A413 airfoil with the 'shape being that predicted by the 
analysis in fiqure 38. The tests were cond ucted in the 6 
by 22 inch transonic wind tunnel located at The Ohio State 
University’s Aeronautical and Astronautical Researoh 
Laboratory. Pour different configurations were tested to 
separate out the roughness and shape effects as are done 
in the analysis. Complete details of the. experiment can 
be found in Appendix A. Here the data are compared, to the 
analysis. 

No detailed pressure data can be found in the 
literature for airfoils with ice shapes, real or 
simulated. Even the most recent, work by the Soviet* 
Swedish group [ 2 ] on simulated ice shapes contains no 
airfoil pressure distributions. These data are necessary 
for a detailed evaluation of the airfoil analysis code. 
Figure 43 shows the measured and predicted Cp distribution 
on the clean airfoil. Here the comparison is made at a 
lift coefficient ol 0.52, Peynolds number based on chord 
l«vigth of thr«‘o million, and larh number of 0.40. The 
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presruTH distribution predicted by the Epplor code is very 
close to that measured in the tunnel. The leading edge 
discontinuity on both the upper and lower surfaces is 
"predicted, although the upper surface is off somewhat in 
magnitude-. The rest of the pressure distribution also 
agrees well. A slight deviation is seen near the trailing 
edge where the boundary layer thicKness affects the 
pressures. This is not a ccou nted for in the current 
version of the Eppler analysis code. 

Figure 44 compares the measured av-d predicted pressaire 
distributions, on the airfoil with the s.^niulated rime ice 
shape. The lift coefficient for this comparison is 0.45. 
The most noticeable feature of the experimental Cp 
distribution are the discontinuous pressure spikes o"n the 
upper and lower surface of the leading edge. These spikes 
are predicted fairly well by the analysis. The presence 
of the spikes will cause early boundary layer transition 
and probably the formation of loaoing edge separation 
bubbles. Therefore the ability of the airfoil code to 
accurately predict this pressure distribution is the first 
stei) toward the accurate analysis of airfoils with rime 
iop . 

The comparison between the measured and predicted lift 
coefficients is shown in figure 45. The predicted angle 
of zero lift compares very well while the lift curve slop? 
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the smooth value of 0.0045. Considering the difficulty of 
the analysis, this represents a reasonabl e comparison. 

Note also tha_t_th£ theory is based on actual iced airfoil 
data and a simulated ice shape was tested. Therefore, the 
error-may be due in part to the way in which the ice shape 
was simulated. This error, in the simulation can not be 
dtiteriuined from these tests, and it suggests that an 
■'experimental prooram is needed to develop ice simulation 
techniques . 

The analytical method for calculating iced airfoil 

performance has been compared to actual airfoil icing 

tests. The predicted ice shape of figure 34 has been 

analyzed and the results are shown in figures 47 and 48. 

The airfoil used is the modified NACA 64-215. The lift 

coefficient curve, figure 47, shows the expected reduction 

in Co due to a leadina edge bubble. Unfortunately no 
'max 

lift coefficient data was tahen on the actual iced airfoil 
to be used for comparison. This reduction in maximum lift 
coefficient does however seem reasonable when compared to 
similar airfoil results. 

The an a iy tica ,lly predicted drag polars for both the 
clean and iced airfoils are shown in figure 4 0. Here 
oxn or i men t a 1 values of T' ag coefficient at 0./ degrees 
angle of attack are available for the clean and iced 
airfoil [68]. Since no ice roughness was reported, the 
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results are shown for values of k/c of 0,001 to 0.005 
which bracket the usual ranqe of rime ice roughness. Here 
the comparison between theory and experiment is very good, 
especially the increnient in i-he. drag due to the ice. 

Again the clean value is calculated using the Eppler 
proaram with his roughness correction and the increase in 
drag is based on eg'. (41) . 

The NACA 65A004 airfoil has been analyzed using the 
rime ice shane predicted in figure 42. The predicted drag 
polar and the measured values are shown in figure 49. 

Here again the experimental drag values are only available 
at one angle of attack. The analysis does an excellent 
job of predicting the drag increase for values in the 
cruise ranac. 

The effect of the ice- Shape on the maximum lift 

coefficient is very unusual for this particular airfoil. 

As seen in figure 42 the ice shape forms a leading edge 

flap for this thin airfoil. The measured increase for 

thirs- case is approximately 23 percent while the analysis 

shows a 12 j'ercent increase in maximum lift -coe-ff-icient . 

Although numerically this comparison may seem less than 

desireable, it actually lends a great deal of confidence 

to the method. The eSAOQh is a very severe test of the 

analysis since the airfoil is so thin. To predict an 

increase in' which is conservative demonstrates that 

'^max 
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the leadina edqe region in being handled correctly by the 
analysis . 

The areodynamic analysis has ~bee-u_xi o . ni . p ared to both 
simulated and actual rime ice on three very different 
airfoil sections. All the_rjasjuIJhSj_ both lilt and drag, 
have compared very well considering the dificulty in 
performing the analysis. The method for the aerodynamic 
analysis of airfoils with rime ice presented here has been, 
shown to be a reliable procedure. Hopefully the empirical 
corrections to the drag predictions can eventually be 
replaced by analytical methods when they become available. 
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VII. SUriMARV AND CONCLUSIONS 

A methodolof^y has been developed to predict the growth 
of rime ice, and the resulting aerodynami.c_4ienalty , oii 
unprotected airfoil surfaces. This, method has for the 
first time included the time effects into the icing 
analysis. A large portion of this study was involved in 
the numerical formulation of the problem for — 

computer solution. Howevex, the derivation of two new 
similarity parameters was primarily an analytical 
exercise, while some experimental work was performed in a 
wind tunnel evaluation of the aerodYnamic analysis. 

The calculation of water droplet trajectories was 
performed by a step inte<jration of the governing stiff 
system of ordinary differential equations. The required 
flowfield was provided by a modified Theodorsen method. 
Although calculations of this type have been performed 
earlier, by usinq state of-the-art computational 
facilities and nUmericai proeeedures a large improvement 
has been made. The present procedure was faster, more 
accurate, and more generally applicable than earlier 


methods . 


I 








An in depth analysis of the governing differential 
equation has lead to a simplified similarity parameter for 
the problem. By using a reduced form of the droplet drag 
equation the two sitnilarity parameters, Ry and K, were 
combined into a single parameter, K, the trajectory 
similarity parameter. This greatly simplified the 
analysis. 

By making a further simplification to the droplet drag 
equation the modifieo inertia parameter, Kq, first 
suggested bv Langmuir, was deriJteji in the same— manner . As 
a result of this analysis a closed form solution was found 
for Kq. This was the first derivation of Kq from the 
governing differential €>quation-, and the first time a 
closed form solution for it has ever been found. 
Experimental and numerical results have been presented in 
support of Kq and K. The new trajectory similarity 
parameter has been found to be superior to Kq, especially 
in scaling applications. 

Using the results of droplet trajectory calculations 
rime ice shapes have been predicted. In the derivation of 
these Gouations a similarity parameter has been 
identified, the accumulation parameter, Ac. For a given 
geometry and K the accumulation parameter governs the 
growth of rime ice on airfoils. 
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As rinip ice builds up on an airfoil leading edge the 
effective airfoil shape becftn).es a ^function of time. This 
then results in the surface flux o f imping ing water 
droplets also being a function of time. The present 
method has included these effects into the ice shape 
prediction. A time -stepping procedure was employed where 
the airfoil geometry, flowfield, and droplet impingement 
efficiencies were updated periodically during the ice 
accretion proc ess . Comi^arison of predicted rime ice 
shapes to those measured in a icina wind tunnel compared 
well. A significant improvement was seen in the 
theoretical shapes when the time-stepping procedure was 
used . 

The t-i me— s tenpi ng procedure has provided insight into 
the ice accretion process. Some researchers have 
suqaested that the ire actually grows out from the surface 
tangent to the incoming droplet trajectories. This 
tanaent ice growth has been shown to be merely an 
aonroximation to the time effects wiijaxfi-JLit g.. growth 

out normal to the surface was used. With the importance 
of including time effects in the rime icing analysis 
demonstrated, the method is expected to provide similar 
improvements to glaze ict^ -predictions. 

The aerodynamic effects of rime ice accretions on 
airfoils includes a reduction in maximum lift coefficient 
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and an increase in draq. Eailit?r methods for predicting 
the degradation in airfoil performance with ice relied 
totally upon empirical correlations. These methods, 
h5wever, dealt only with the changes in drag and were 
based on initial icing rates. The present method for 
evaluating the iced airfoil performance was based on an 
analytical analysis of the resulting airfoil shape after 
ice accretion. The mehtod postulated that the aerodynamic 
effects of riiue ice were due to; 1) the surface roughness 
of the ice, and 2) the change in leading edge geometry due 
to the smooth ice shape. These two mechanisms were then 
handled separately by the analysis. 

The smootfi ice shape was analyzed using existing 
airfoil analysis codes. The surface roughness effect was 
handled by correcting the analytical results based on an 
empirical equation which was developed here. 

Since no detailed aerodynamic data on an airfoil with 
rime ice was available, wind tunnel tests on an airfoil 
wi th simulated ri me i re were cond iictfid-.^ The experiment 
identified the effect surface roughness and ice shape have 
on airfoil performance. In addition to lift and drag 
data, these tests generated the first detailed pressure 
measurements ever taken on an airfoil with simulated ice. 
The predicted pressure distributions compared well with 
t-he experimental— results as did the values for and Cd, 
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The aerodynamic analysis was verified further using values 
of lift and dran fr-oja icing wind tunnel tests of actual 
ice accretions. 

The present study has identified areas where 
additional work is noedeit. The analytical method could be 
improved by either removing the need to smooth the shape 
or iinjiroving the smoothing procedure. This would increase 
the accuracy of the. ice shape prediction and allow smaller 
step sizes. In addition better information on the ice 
density wovilO greatly improve the me44^ed . Future 
analytical research on rough airfoil drag could remove the 
need to use an empirical" drag correlation. Experimentally 
the need is to expand the old, and very limited, data base 
in terns of accurate ice shapes, ice densities, and 
airfoil aerodynamic performance penalties. However the 
most serioiis need is to extend this work to the glaze ice 
case where a flowfie-ld with large zones of separated flow 
must be accurately predicted. 

In summary the rime icing— methodology prese_nte^ here 
has advanced the state-of-the-art in four major areas. 
First, the effects of time on the ice accretion process 
have been i-ncludeo in the analysis. By using the tlme- 
stepning method- very accurate rime ice shapes can be 
predicted. Secoiul, an aerodynamic analysis has been 
formulated which is based on th<* actual iced airfoil 
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qwau’try. On like early methods which estimated Cd from 
only initial icina rates, thi-s mothou predicts and Cd 
from the new airfoil geometry with some ettpirical 
corrections . . . .. 

The third major contribution came from the wind tunnel 
test of the simulates ice shape. Here for the first time 
detailed aerodviiaraic data, including surface pressure 
distributions., were taken on an„airfoil with simulated 
rime ice. The data provided a great deal of insight into 
the problem and an excellent test case for the present, and 
tor future aerodynamic analysis. The similarity analysis 
provided the final contribution. Two new parameters, K., 
the trajectory similarity parameter, and, Ac, the 
accumulation parameter have been derived and shown to 
govern the accretion of rime ice on airfoils. In 
addition, K^, the modified inertia parameter has been 
derived from the governing differential equat-ion and the 
first closed form solution for K(j_has been presented. 
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FIGURE 10. FLOWCHART FOR THE FaME’ ICE METHODOLOGY 
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FIGURE 1 2. DROPLET INITIAL Y COORDINATE IN THE FREE STREAM AS A 
FUNCTION. OF IMPINGEMENT POINT ON THE AIRFOIL 
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FIGURE 17. EMPIRICAL FIT TO ICED AIRFOIL DATA 
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FIGURE 18 COMPARISON OF THE PRESENT METHOD TO NACA TN 3586 FOR 
IMPINGEMENT EFFICIENCY ON A NACA 65A004 AIRFOIL 
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FIGURE 19. COMPARISON OF THE PRESENT METHOD TO THAT OF LOZOWSKI FOR THE' I 
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FIGURE 20- PRESENT METHOD FOR CALCULATING IMPINGEMENT EFFICIENCY COMPARED 
TO LOZOWSKI'S WITM DRAG DUE TO DROPLET ACCELERATION 
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FIGURE 22. EFFECT OF DROPLET DIAMETER ON IMPINGEMENT EFFICIENCY 
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FIGURE 30 


DROPLET TRAJECTORIES ABOUT A NACA 0012 AIRFOIL AT ZERO AND FIVE 
DEGREES ANGLE OF ATTACK 
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FIGURE 34. PREDICTED RIME ICE SHAPES FOR- THREE 
FIVE MINUTE TIME STEPS 
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FIGTOE 35. lOTINGEMEOT EFFICIENCY CURVES FOR THE MODIFIED NACA 64-215 AIRFOIL 
ANALYZED USING THREE FIVE MINUTE TIME STEPS 
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FIGURE 37 IIIPINGEMENT EFFICIENCY CURVES FOR JTHE MODIFIED NACA 54-215 
AIRFOIL ANALYZED USING SIX 2.5 MINUTE TIME STEPS i 
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FIGURE 38. 


EFFECT OF THE NUMBER OF TIME STEPS ON THE 
PREDICTED RIME ICE SHAPE FOR A NACA 65A413 
AIRFOIL 




FIGURE 39- 


IIIPINGEMENT EFFICIENCY CURVES FOR THE TWO 
TIME STEP ICE CALCULATION 
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FIGURE 40- IMPINGEMENT EFFICIENCY CURVES FOR THE FOUR TIME STEP 
ICE SHAPE CALCULATION 
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FIGURE 41. THEORETiCAI. RIME ICE SHAPE COMPARED TO 
EXPERIMENT FOR THE MODIFIED NACA 64-215 
AIRFOIL 
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— - — 1 5 MIN STEP, E * .153 




-i . 


THEORY 



IGURE 43. THEORETICAL AND EXPERIJIENTAL PRESSURE DISTRIBUTION 
FOR THE CLEAN NACA 65A413 AIRFOIL 
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FIGURE 44. THEORETICAL AND EXPERIMENTAL PRESSURE 

DISTRIBUTION FOR THE NACA 65A413 AIRFOIL 
WITH SIMULATED RIME ICE 
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FIGURE 48. EXPERIMENTAL AND THEORETICAL DRAG POLARS 

FOR THE MODIFIED NACA 64-215 AIRFOIL CLEAN 
AND WITH ICE 



155 


rppfndix 


ha- been conducted at The. Ohi<^ 

,, experimental program ha. be. 

. „itVs Aeronautical and Astronautxca.l , 

Characteristic* ci 

Laboratory to determin tunnel test 

r^atL 

,,,, to ha ara^ „„i, to. 

la'te. rlae ic.aeroa.na.ic lata, hat also to 
,„orato ,„al,tical .eihod that the 

test the h.ia,thests asei 
.„ects of ice shape a„a roa,hness can 

.... ■•“• ■••■ ■“•' 

confianrationn: 

1.) Clean airfoil (baselxn^)^^ rouqhnesB 

liffSl :lu ooooth ri.e/- -rreaiin, ei,e 

“•> *“arhes"/Sara ,si.ala«^ ri.e ice, 

^-haracteristics of each 
♦ Ko aprodYnamic characrei- 

the ettects ot sar.ace roao.hhes* ani ice 

cohfiaarattoh. the 

Shape can be determine . tioure IS can be 

,.ecK on the ACQ prediction of 
base-lin-e a Brumby, 

4^ as weXl as a checX QD_the ACg,^ax 
tnaae as 


156 


figure '(3. Tho> results of models 3 and 4 compared to the 

baseline will provide verification of the Co analysis. 

■^max 

Finally the tests of model 4 will verify the entire 
theo retica 1 method. 

Experimental Facility 

The experimental facility used in this study was the 
OSD 6 by 22 Transonic Ai-rfoil wind Tunnel [69], The 
tunnel is designed for two dimensional testing with a test 
section six inches wide, twenty-two inches high, and 
forty-four inches long.- The side walls are solid, while 
the top and bottom walls of the tunnel are perforated with 
a porosity of ten percent. The tunnel operates in a 
blowdown mode with the F!ach number controlled by a choRo 
downstream of the test section. Mach numbers from 0.2 to 
1.1 arc available. The total pressure in the stagnation 
chamber is varied to control the Reynolds number and 
provide a range of 1.5 to 33 million per foot. 

Lift and moment coefficient data are normally taRen 
usino model static pressure taps. Pressure measurements 
are made with a Scanivalve, trapped volume system which is 
sampled with a transducer after the tunnel is shut down. 
Drag data are taken using a wake survey probe which 
traverses the wake recording the staonation pressure 
deficit. The data collected is digiti 2 ed and stored on 
n>\onetic tape in the Harris SLASH 6 Digital Computational 


157 


Facility [70] of the Aeronautical and Astronaut ical 
Rese-arch Laboratory. The data is then reduct?d to 
coefficient form [71] and output as quick look data on a 
CRT display or hard copy printed and plotted. 

The interference effects in the OSO 6 by 22 Hind 
Tunnel have been investigated [69], Confinement 
interference, spanwise interference from the side walls, 
and flow quality have been evaluated. .The correction 
required for six inch chord models-has been shown to be 
negligible. The correction to tne angle of attack is on 
the order of 0 . 17- degrees per unit Cj^. Since this test 
will use a six inch chord airfoil, no corrections need be 
made to the data. 

Airfoil Model 

A NACA 65A413 airfoil section was selected for the 
experiment. The model usod in the wind tunnel was a brass 
model of six inch chord and six inch span, figure A-1. 

The original airfoil model was instrumented with 96 static 
pressure taps of which 42 were used in the data reduction. 
The trailing edge tap is located -on the sidewall due to 
the physical constraints. 

The rinte shape which was simulated was that predicted 
by the time-stepping analysis of figure 37. A comparison 
of the predicted shape and the shape used on the tunnel 
model is given in figure A-2. Note that~since the 
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objective of the test was to generate baseline data to 
validate the analysis, the accurate reproduction of the 
predicted shape is not required. All that is required is 
that the ice shape siaulated be a representative geometry 
and that it be adequately documented for the analytical 
com parisen . 

A schematic of the airfoil model with the simulated 
rime ice shape is shown in figures A-3 and A-4. The rime 
ice shape was simulated by adding a 0.145 inch outside 
diameter tube to the airfoil leading edge. The mounting 
blocks were drilled to allow the tube to extend out of the 
tunnel on both sides. The center section of the tube was 
rPDlaced by a solid rod which was drilled through to pick 
up the existing leading edge pressure tap. The tube, now 
plugged in the center, was used to add an additional tap 
on the leading edge uppor annd lower surfaces, figure A-4. 
The ice shape was completed by building up the area 
between the tube and the airfoil until the desired shape 
was reached. Care was taken to ensure that the affected 
airfoil pressure taps were exteijded up through this xeq ion 
to the new airfoil surface, A photograph of the airfoil 
mcKlel with the simulated rime ice shape is shown in figure 
A-5. 

Roughness was added to the moael for configurations 2 
a-n4-4. This roviqhness was intended to simulate the actual 
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rouahness on a rime ice shape. Rime ice surface roughness 
is typically in the range of k/c » 0.001 to 0.005. 
Carborundum grit with an average size of 0.015 inches was 
used. This scales to a k/c of 0.0025 for a six inch chord 
model . 

The grit was applied by first coating the surface with 
Krylon clear acrylic spray to provide the adloes.l«.e.. The 
roughness— elements were then applied to the surface and 
two or three coats of acrylic were applied to ensure that 
the particles were firmly adhered to the surface. The 
rouqhness was applied ts the leading edge of the airfoil 
upper and lower surface back to three percent airfoil 
chord for both configurations 2 and u. The roughness 
elements were distributed randomly at a concentration of 
about 250 per square inch of surface area, figure ft-6. 
Results and Piscussioii 

The airfoil section selected is typical of that 
curren-tly in use on general aviation and business 
aircraft. To simulate actual operating condi tions, .a Mach 
number of 0.40 and a Reynolds number based on chord length 
of 3 million were chosen for the cruise case. These 
conditions were used in testing the airfoil at angles of 
attack of eight degrees and less. To determine the 
maximum lift coefficient, conditions more typical of a 
landing approach wer.e. used. For angles, of attack greater 
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than eight degrees a flach number of 0.23 and Reynolds 

numner of 2 million were used, _ 

Pressure distributions for the clean airfoil and with 
surface roughness added, configurations 1 ami 2, are shown 
in figure A-7. Here both air.foils are at two degrees 
angle of attach. Both curves are similar, however the 
areas, which give the model lift coefficient, are 
different. The model with roughness e-X^^eriences a 
decrease in lift over the clean model. This is pxobably 
due to the effect of the roughness on the boundary layer. 
The roughness results in a thicXer boundary layer at the 
trailing edge upper surface and therefore a larger 
displacement thickness. This effectively removes camber 
from the airfoil and decreases the lift, shifting to a 

more positive value. 

Note also tlve reduced pressure recovery at the 
trailing edge for the rough airfoil. This suggests 
increased drag which cauv be easily scfcn in the wake 
deficit plots of figure A-8. Here the roughened airfoil 
has a larger velocity deficit, and therefore more drag. 
This result is in good agreement with earlier experimental 
work showing increased drag with surface roughness. 

The airfoil with simulated rime ice experiences 
somewhat different surface pnessures near the leading 
edge, figure R-9. Note that here the chord lengJLh is, _ 
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based on the iced airfoil chord, 1.024 times the original 
chord. The alt portion of the Cp distribution is similar 
to the no ice case, however the leading edge region is 
altered by the ice shape. Pressure spikes, severe 
discontinuities in the pressure distribution, occur on the 
upper and -lower surfaces where the ice shape joins the 
airfoil contour. For this ice shape this reprei^njts a 
discontinuity in the second derivative of the surface 
shape. These spikes wer*; detected by the two additional 
pressure taps installed in the tube which forms the 
•leading edge shape of the simulated ice. This 
demonstrates the importance of the installation of. 
pressure taps in simulated ice shapes. The effects of the 
pressure spikes will be seen to be more serious at higher 
angles of attack. 

Figure A-10 shows the lift coefficient as a function 


of angle of attack for all four airfoil configurations. 
Configurations 2 through 4 all have approximately the same 
effect on the lift coefficient. Those changes are a shift 
in ttLO ^ sizeable decrease in good 

agreement between the fbe smooth and rough rime 

ice shape suggests that the stall is caused by the shape. 
As seen in figure A--^, the severe pressure gradients near 
the leading edge probably lead to a leading edge 
separation bubble. This accounts for the early separation 
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at hiqher angles of attack. 

The reduced for the airfoil with only surface 

^max 

roughness is due to a different mechanisro. Here the 

roughness cau-ses a thickening of the boundary layer and 

decreases pressure recovery at the trailing edge. This 

ultiirately leads to early trailing edge separation which. 

moves forward as a increses to cause the reduction in 

maximum lift coefficient. The apparent agreement in Cp 

'^max 

for configurations 2 and 4 is due to the particular k/c 
and rime ice shape chosen, and should not be interpreted 
as a general trend. 

Figure A-11 shows the measured drag polars for all 
four conf iauration.s . The smooth airfoil is seen to have a 
minimum drag coeficient of about 0.086. This is well 
above the laminar "drag bucket” values expected for an 
airfoil of this type. The brass airfoil model used was 
slightly tarnished and therefore did not have the surface 
finish necessary to permit long lamiaax. runs- 

An increment in dra-g was seen due to the addition of 
surface roughness. This drag increase was certainly 
expected and is of a reasonable magnitude. The reason for 
the apparent agreement between the rough airfoil and the 
smooth ice shape, configurations 2 and 3, is not obvious. 
Most- likely this is not a general result, but again merely 
a coincidence resulting from the roughness and the rime 
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ice geometry chosen. 

An additional draq increment was measured when 
roughness elements were added to the rime ice shape, 
configuration 4. This, is the simulated rime ice shape. 

This increase in drag contrasts the maximum lift 
coefficient case were configurations 3 and 4 behaved 
similarly. Theretore , "while ice shape alone is sufficient 
to determine » the surface roughness of the ice nurst_ 

be modeled to simulate accurately the total drag increase 
due to airfoil rime icing. 

The moment coefficient about the quarter point of the 
original airfoil is plotted as a function o, lift 
coefficient in figure A -12. Configurations 2 through 4 
all show a reduction in the nose down pitching moment when 
compared to the Glean model. The leading edge roughness 
as described before thickens the boundary layer and 
unloads the aft portion of the airfoil section. Therefore 
the nose do-wn pitching moment is reduced. The smooth ice 
shape adds area in front of the nose providlag_BUirsJ_lMaS£ 
up moment, explainina conf iauration three's reduction in 
nose down pitching moment. The rough ice shape combines 
the two above effects, resulting in a slightly larger 
rfjduction in nose down pitching moment than that 
experienced by the shape alone. 
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This not only provided data tc^ verify the 

analysis, but has demonstrated the feasability of 
performing siiaulated iced airfoil tests in a small scale 
wind tunnel facility. The data show the expected results 
of- decreased maximum lift coefficient and increased drag 
with the simulated ice shape. Jn addition a reduction in 
nose down pitching moment was measured with simulat_e d r ime 
ice. The pressure distribution: measured for the airfoil 
with simulated rime ice are believed to be the first such 
data published. These Cj:> plots provide insight into the 
physicjil phenomena and detailed information to be used to 
evaluate and refine current analytical methods. 
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FIGURE A- 3- TOP 


VIEW OF SIMULATED RIME ICE AIRFOIL MODEL 
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FIGURE A- 4. 


CROSS SECTION OF SIMULATED RIME ICE AIRFOIL MODEL 
SHOWING PRESSURE TAP LOCATIONS 




FIQURE A-5. WIND TUNNEL MODEL OF AN NACA 65A413 
AIRFOIL WITH SIMULATED RIME ICE 
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FIGURE A- 6. 


LEADING EDGE OF SIMULATED RIME ICE MODEL 
SHOWING SURFACE ROUGHNESS 
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VELOCITY RATIO, u/U 


FIGURE A- 8. 


VELOCITY RATIO IN THE WAKE OF THE NACA ^5A413 AIRFOIL 
CLEAN AND WITH LEADING EDGE ROUGHNESS 




EOTERIMENTAL PRESSURE DISTRIBUTION OF AN NACA 
ccA/.n ATurnTT UTTH qiMULATED RIME ICE 


FIGURE A- 9. 
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FIGURE A-10. LIFT COEFFICIENT AS A FUNCTION OF ANGLE 
OF ATTACK FOR THE NACA 65A413 AIRFOIL 
WITH LEADING EDGE MODIFICATIONS 
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FIGURE A- 11. 


DRAG POLARS FROM THE WIND TUNNEL TESTS OF THE 
NACA 65A413 AIRFOIL WITH LEADING EDGE MODIFICA 
TIONS 



FIGURE A- 12. 


MOMENT COEFFICIENT AS 
THE WIND TUNNEL TESTS 
EDGE MODIFICATIONS 
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